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The g-model, a random walk model rich in behaviour and applications, is investigated. We 
introduce and motivate the g-model via its application proposed by Coppersmith et al. to the flow 
of stress through granular matter at rest. For a special value of its parameters the q- model has a 
critical point that we analyse. To characterise the critical point we imagine that a uniform load has 
been applied to the top of the granular medium and we study the evolution with depth of fluctuations 
in the distribution of load. Close to the critical point explicit calculation reveals that the evolution 
of load exhibits scaling behaviour analogous to thermodynamic critical phenomena. The critical 
behaviour is remarkably tractable: the harvest of analytic results includes scaling functions that 
describe the evolution of the variance of the load distribution close to the critical point and of the 
entire load distribution right at the critical point, values of the associated critical exponents, and 
determination of the upper critical dimension. These results are of intrinsic interest as a tractable 
example of a random critical point. Of the many applications of the q-model, the critical behaviour 
is particularly relevant to network models of river basins, as we briefly discuss. Finally we discuss 
circumstances under which quantum network models that describe the surface electronic states of a 
quantum Hall multilayer can be mapped onto the classical g-model. For mesoscopic multilayers of 
finite circumference the mapping fails; instead a mapping to a ferromagnetic supersymmetric spin 
chain has proved fruitful. We discuss aspects of the superspin mapping and give a new elementary 
derivation of it making use of operator rather than functional methods. 

PACS: 



I. INTRODUCTION 

It is fortunate that in physics the same equations some- 
times arise in contexts that are apparently very different. 
Feynman illustrates this through elementary examples in 
his introductory lectures on physics to impart the les- 
son that the "same equations have the same solutions" 
0. Our purpose is to study a model, recently dubbed 
the q— model, that provides another such instance. The 
q— model has been used to describe the merging of trib- 
utaries to form rivers the aggregation of diffusing 
charges ji); the flow of stress in a granular medium jjj; 
and can be mapped onto the abelian sandpile, a model 
studied in context of self-organised criticality ||. It is 
also closely related to models that describe the surface of 
a quantum Hall multilayer and passive scalar turbu- 
lence UK . Here we focus on the application to granular 
matter, river networks and the quantum Hall multilayer. 

Granular matter exhibits fascinating behaviour that 
is little understood jlQ] . Examples of granular matter 
include sand, powders and agricultural grains stored in 
silos. An important problem is the propagation of stress 
through a granular medium at rest. This has been stud- 
ied by ingenious experiments, in which a vertical load is 
applied to an amorphous pack of beads, and the loads on 
the beads in the top and bottom layers are recorded using 
carbon paper ]ll| , p^ |. Such experiments yield the distri- 
bution of load on the beads and reveal that there are no 
horizontal correlations in load even amongst neighbour- 



ing beads. The g-model was introduced by Coppersmith 
and coworkers to account for the distribution of load [EJ . 
As we shall see, it also correctly predicts the lack of hor- 
izontal correlation. 

For simplicity we describe the g-model in a plane. 
Since the vertical and horizontal directions are treated 
asymmetrically we call this the 1+1 dimensional g-model. 
The extension to 2+1 dimensions (relevant to experi- 
ments on bead packs) and higher, is straightforward and 
is discussed in section V. In the g-model it is assumed 
that the beads sit on a regular lattice shown in fig 1. 
The location of the beads is specified by the co-ordinates 
t (the depth of the layer) and n (the location of the bead 
within the layer). Note that n takes only even values for 
t even; only odd, for t odd. Each bead is assumed to 
be supported by its two nearest neighbours in the layer 
directly below. More precisely, it is assumed that a ran- 
dom fraction /„ (t) of the load of bead (n, t) is supported 
by the neighbour to the left, bead (n — 1, t + 1); the re- 
mainder, 1 — f n (t), by the neighbour to the right, bead 
(n + l,t + 1). Denoting the load on a bead w and its 
weight I we may write 

W n {t) = Wn-l(t - 1)[1 - /„_i(i - 1)] 

+w n+1 (t-l)f n+1 (t-l)+I n (t). (1) 

The content of eq (1) is that the load on a bead is the 
sum of the loads transmitted to it by its neighbours in 
the layer above plus its own weight. The last term in eq 
(1) is called the injection term. Once the fractions are 
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specified, a given load on the top layer can be propagated 
downward by use of eq (1). 

In the g-model it is assumed that the fractions are in- 
dependent, identically distributed random variables. The 
distribution is assumed to be symmetric about / = 1/2 to 
avoid introducing a horizontal drift to the flow of stress; 
in other words it is assumed P(f) = P(l — /). There is 
no other restriction. Thus the g-models really constitute 
an enormous family of models corresponding to different 
symmetric distributions P(f). To fully specify a particu- 
lar model it is necessary to choose the distribution P(f). 
One obvious possibility is to take P(f) to be uniformly 
distributed over the unit interval; another is to assume 
that the fractions must be or 1 with equal probability. 
The latter is called the singular distribution. 

Mathematically, the g-model is a problem of random 
walkers that coalesce upon contact and fission sponta- 
neously. The singular distribution corresponds to the 
case that the walkers coalesce but do not fission. 

Coppersmith et at argued that, neglecting injection, 
at sufficient depth the distribution of load would attain a 
steady state Q. They studied U(w,t — > oo), the proba- 
bility distribution of load on beads in a sufficiently deep 
layer. For almost all distributions P(f), except the singu- 
lar distribution, they concluded that U(w, t — * oo) decays 
exponentially for large w. This agrees with experiment 
and constitutes an important success of the g-model. For 
the singular distribution, Coppersmith et al. argued that 
H(w, oo) follows a power law. Hence they conjectured 
that the singular distribution constitutes a critical point 
in the family of g-models. A major goal of this paper is to 
make this analogy to thermodynamic critical phenomena 
precise by detailed analysis of the critical point. 

In spite of the success mentioned above the g-model 
cannot be considered a complete theory of stress propa- 
gation in granular matter. This is clear both empirically 
and on grounds of internal consistency. Since the publica- 
tion of the g-model, interesting new ideas on the subject 
of stress flow have appeared ||[l3| ^5[ , but in this paper 
we restrict attention to the g-modcl. This seems justi- 
fied because the g-model does capture some elements of 
the physics correctly and because it exhibits non-trivial 
critical behaviour that is interesting in its own right. 



Further motivation to study the g-model and particu- 
larly its critical point comes from hydrology. To make 
contact with that subject consider a singular g-model 
with zero injection and imagine that only a few beads in 
the top layer are loaded. The load then zig-zags down- 
wards, perhaps along the lines shown in fig 2. If we in- 
terpret these lines as tributaries merging to form a river 
we arrive at Scheidegger's model Q which appeared in 
the hydrology literature more than thirty years ago . 
Networks of tributaries in river basins are known empir- 
ically to be scale invariant structures that obey a variety 
of power laws. Scheidegger networks too obey these laws 
and are in this statistical sense extremely realistic repre- 
sentations of river basins. An excellent discussion of river 
basin power laws is given in refs E8 19 1 . Ref 1 18 1 presents 



some discussion of data; ref |r| provides a detailed com- 
parison between real and Scheidegger networks. 

Here we wish to point out that non-singular g-models 
too can be interpreted as models of river networks. For 
example, consider a model in which the fractions can take 
only the values 0, 1/2 and 1 with probability (1 — S)/2, 5 
and (1 — S)/2 respectively. This model reduces to Schei- 
degger's as 8 — > 0. It produces networks similar to Schei- 
degger's except that occasionally streams split to form 
distributaries. Thus this network is topologically dis- 
tinct from Scheidegger networks. More significantly, as 
we show below, a network with non-zero S is not scale- 
invariant. This is reminiscent of a river network model 
studied by Narayan and Fisher po| . In their "rocky- 
river" model too the network is not scale invariant ex- 
cept if a model parameter is tuned to a special (criti- 
cal) value. Effectively this tuning parameter also con- 
trols river splitting. Taken together, these results suggest 
that river splitting is a relevant perturbation that spoils 
the scale invariant structure of networks. In this paper 
we concentrate on showing that g-model networks with 
river splitting are not scale invariant. We do not explore 
whether such non-scale invariant networks are realised 
in nature (for further discussion and speculation in this 
regard, however, see section VII). 

A quantum Hall multilayer consists of layers of two- 
dimensional electron gases stacked vertically. Multilay- 
ers can be realised by fabricating an appropriate GaAs 



1 Parenthetically we note that Scheidegger's model is purely 
descriptive in the sense that it is a recipe to draw statis- 
tically realistic networks. Somewhat different in spirit are 
models that seek to represent physical processes, sometimes 
very crudely, by which the network forms. Two examples of 
such models in the recent Physics literature are refs 
The model of Leheny and Nagel for example describes an 
apocalyptic lattice world with discrete time. Each time step 
brings precipitation, and in its wake, erosion and avalanches. 
Realistic networks result. 
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heterostructure fl2l|. They are also realised naturally in 
some organic salts. In a quantum Hall multilayer a suf- 
ficiently large magnetic field is applied perpendicular to 
the layers so that the lowest Landau level in each layer 
is fully occupied. Under this circumstance the only im- 
portant electronic states in each layer are the chiral edge 
states that propagate in one direction only as shown in 
Fig 3(a). These edge states are coupled by tunneling be- 
tween layers. Thus the surface of a multilayer is covered 
by a chiral sheath of coupled edge states. These surface 
states control the electrical transport properties of the 
multilayer. A central question from a quantum transport 
point of view is whether these surface states are localized 
or extended in the direction of the field [[^, ^2|j23]| . 

Fig 3(b) shows a network model of the multilayer sur- 
face introduced by Saul, Kardar and Read Q and stud- 
ied by many authors subsequently. In this model it is 
assumed that tunneling between edges takes place only 
at discrete nodes (dashed vertical lines in Fig 3b) that 
appear at regular intervals along an edge. The edges are 
separated by nodes into horizontal segments called links. 
The wavefunction has a definite value on each link. Each 
node is visited by two incoming links and by two outgo- 
ing links. Each node is characterised by a 2 x 2 S-matrix 
that relates the wavefunction on the outgoing links to 
the incoming amplitudes. Once the S- matrices are spec- 
ified, given the wavefunction through a vertical slice, we 
can propagate it to the right. The S- matrices are chosen 
at random from some ensemble to incorporate the effect 
of disorder. To fully specify the model it is necessary to 
choose a distribution for the S-matrices. Periodic bound- 
ary conditions are imposed in the horizontal direction [Q . 

The directed network model above is quantum mechan- 
ical but in the limit of infinite circumference and for a 
special choice of disorder, Saul, Kardar and Read have 
shown that it reduces to a classical model, the g-model 
with uniform distribution of fractions and zero injection 
. In section VI we discuss some respects in which more 
generic models of the multilayer surface, that do not re- 
duce to classical models, still do show behaviour similar 
to the g-model |24|,|25|]. At the same time we show that 
in case of finite circumference quantum interference ef- 
fects become important and there is little to be learnt 
from the study of the classical g-model. Instead a map- 
ping to a ferromagnetic supersymmetric spin-chain has 
proved fruitful in this case [ p3"p6j ]. In section VI we dis- 
cuss aspects of this mapping. 

A detailed summary of our results is given in section 
VII. The reader interested in first obtaining an overview 
of the paper or interested only in the results should pro- 
ceed directly to section VII. 



II. CRITICAL BEHAVIOUR IN 1+1 
DIMENSIONS 

Coppersmith et al. analysed the distribution of load 
H(w, t — > oo) at very large depth where presumably a 
steady state is achieved H . Here we study how the dis- 
tribution evolves as a function of depth to this asymptotic 
steady state. We assume that a uniform load is applied 
to the top layer, 

w n (t) = 1 for all n. (2) 

In this section we neglect the weight of each bead (the 
injection term). In partial support of this neglect we note 
that in the experiment of ref jL2| typically a total load of 
7600 N was applied to the bead pack. In comparison we 
estimate that the weight of a single bead was less than a 
mN; of the entire pack, less than 100 N. However, right 
at the critical point injection is a relevant perturbation, 
and at sufficiently large depth must be taken into account 
even if the weight of a single bead is small. We return to 
the effects of injection in section IV. 

To make the problem tractable we study not the en- 
tire distribution H(w,t) but only its lowest non-trivial 
moment. With the neglect of injection it follows that the 
total load on every layer is the same; the g-model dy- 
namics (eq 1) just shuffles this load. Hence the average 
load in layer t 

/>oo 

(w(t)) = / dwwU(w,t) 
Jo 

= 1- (3) 
The lowest non-trivial moment is therefore the variance 

(6w 2 {t))= dww 2 n(w,t) - 1. (4) 
Jo 

Since a uniform load is applied to the top layer the vari- 
ance in that layer vanishes. As the load propagates down- 
ward, the fluctuations must grow and saturate. Our pur- 
pose is to analyse this evolution for different distributions 
P(f), particularly those that are close to the singular dis- 
tribution. 

Right at the critical point the asymptotic distribution 
Tl(w, oo) is believed to be a power law. If we assume that 
it does not have a well defined variance, then by analogy 
to critical phenomena we surmise that close to the critical 
point the variance must diverge as 

(to 2 (t->oo))~^. (5) 

Here d measures the distance of a distribution P(f) from 
the singular distribution; <5 will be defined precisely be- 
low. We also expect that the depth-scale £ CO rr at which 
the steady state is attained will diverge as the critical 
point is approached. Thus 
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£ CO rr is a vertical correlation length that diverges as the 
critical point is approached. Combining eqs (5) and (6) 
we expect that close to the critical point the fluctuations 
must have a scaling form 



(7) 



To be consistent with cq (5) we expect that the scaling 
function !F{u) — > const as u — ► oo. For short times we 
expect that the system should behave as it would at the 
critical point. The 8 dependence should cancel and so we 
expect T(u) ~ u 9 ^ for u < 1 so that (Sw 2 (t)} ~ t 8 ^ at 
the critical point. 

In the remainder of this section we will confirm that 
eq (7) and these inferences are valid. We will determine 
the exponents 6 and tp and the scaling function T{u). 

As an aside to experts we note that it may have been 
more natural to name the exponents 6 — > (3 — t)/ct and 
ip — > vz. These names follow from a more general scal- 
ing hypothesis for the entire distribution (eq 174). How- 
ever in this section we have elected to make the more 
restricted hypothesis eq (7) and to give the exponents 
single letter names taking care to avoid common expo- 
nent names such as a, (3 and v. 



A. Disorder Average 

Consider the correlation function 



(8) 



We assume there are N beads in each layer and we impose 
periodic boundary conditions in the horizontal direction. 
Ultimately we are interested in taking N — > oo. Note 
that m is even for both t = even and t = odd. In terms 
of the correlation function the variance is given by 



(5w 2 (t)) = c (t) - 1. 



(9) 



The correlation function obeys a remarkably simple evo- 
lution equation. This equation can be solved by straight- 
forward classical analysis to yield the evolution of the 
variance. It is not difficult to obtain the entire correla- 
tion function by this method, and thereby obtain infor- 
mation on the horizontal correlation length, but we do 
not attempt this here. 

To analyse the evolution of the correlation function we 
write 

Cm(t + 1) = y^(w n (t + l)w n+m (t + 1)) 

n 

= ^2(Wn+l(t)w n+m+1 (t)f n+1 (t)f n+m+1 (t)) 

n 

+ others. 



^ (Wn+l (t)w n+m+1 (t)) (/„+! {t)f n+m+1 (*)) 

(10) 



N 

n 

+ others. 



To obtain the second line of eq (10) we have used eq (1). 
Four terms result; we have written only one for illustra- 
tion. To obtain the third line it is crucial to observe that 
w n (t) depends only on fractions in the layers above. It 
is not correlated with the fractions in layer t, allowing us 
to factorise the average as shown. 

To perform the average we need information about the 
distribution P(f). By symmetry for any choice of distri- 
bution 



= j\ffP{f) = \. (11) 



For the variance we write 



( ) = 



(12) 



e is a parameter that characterises the distribution P(f). 
For example, e = 1/3 for the uniform distribution. For 
the singular distribution the parameter takes its maxi- 
mum possible value e = 1. Since the fractions for differ- 
ent beads are assumed to be independently distributed 
we conclude 



(/n(*l)/m(*2)) = ^ + jSn^mSt! 



t 2 - 



(13) 



Substituting eq (13) in eq (10) we obtain 

C m (t + 1) = Q + ^m,0^ C m (t) 



others 

2 + !<*m,0 ) C m (t) 



- [ j - |<W I <■„, , 2 (h 

- ( \ - ) Cm I _'(/)■ 



(14) 



In the second line of eq (14) the other terms have been 
unveiled. Recall that m takes even integer values. It is 
convenient to replace m — > m/2 to obtain 

C m (t + 1) = Q + |^m,0^ C m (t) 

^-^m,-ijc m+ i(«). (15) 

Eq (15) is the main result of this subsection. It gov- 
erns the evolution of the correlation function. We wish 
to solve it subject to the initial condition 
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c m (t — > 0) = 1 for all m. 



(16) 



The initial condition follows from the definition of c m (eq 
8) and the assumed uniform load on the top layer. Note 
that the distribution P(f) enters the evolution equation 
only through the parameter e. Since the parameter takes 
its maximum value e = 1 for the singular distribution we 
may define 



S = 1 



(17) 



as the distance of a distribution P(f) from the critical 
point. 



A complication we must negotiate is that H is non- 
Hermitian. According to the standard theory of biorthog- 
onal expansion (briefly recounted in Appendix A) to ex- 
ecute the plan above we must prove that the eigenvec- 
tors of H span the vector space. Then we must find 
the eigenvectors of , called the left eigenvectors of H 
in this context. The eigenvalues of ijt are the complex 
conjugate of the eigenvalues of H . Thus 



(24) 



ip* denotes the left eigenvector with eigenvalue A* . Hav- 
ing completed these tasks we may write the completeness 
relation 



B. Scattering solution 

It is easy to verify that a steady state solution to eq 
(15) is 



co = -, ; 

1 — e 

c n = 1 for n 7^ 0. 



(18) 



Assuming this is the unique steady state towards which 
our initial condition evolves, eq (18) reveals that the 
variance does diverge as the singular distribution is ap- 
proached. Using eq (9) 



(S 2 w(t -> oo)) = 



1 -e 
1 



as e — ► 1. 



(19) 



Comparing eq (5) we see that the exponent 9 = 1. Eq 
(18) also reveals that in steady state the fluctuations 
in load are uncorrelated for all pairs of beads including 
neighbours. This is in agreement with experiment Jl2[ . 

A full solution of evolution dynamics needs more work. 
Schematically eq (15) states 



c(t + 1) = Hc(t). 



(20) 



The strategy we adopt here is to seek the eigenvectors of 



H<P X = \<P X 



(21) 



and to expand the initial correlation vector c(0) in terms 
of the eigenvectors, 



C (0) = ^a A ^ A . 

A 

The correlation vector at depth t is then 



(22) 



(23) 



A 



(25) 



Using eq (25) we conclude that the expansion coefficients 
in eq (22) are determined by the left eigenvectors: 



ax 



(26) 



Implementing the plan we first write the eigenvalue 
equation for H 



2^ 

0V 



1 

"4" 

1 

-<t 

4 ' 



w r+l 



'-2 



1 

f 4< 
1 - e 



V-i 



Mr for \r\ > 2; 



-<f) Q = X^_i, 



1 



+ + \& = Hi- 



(27) 



Note that for e = eq (27) may be interpreted as the 
Schrodinger equation for a free particle on a tightbinding 
lattice, familiar from elementary solid state physics. For 
non-zero e the particle may be viewed as scattering off 
a (non-Hermitian) barrier at the origin. Thus we seek a 
solution of the scattering form 



0(+)fc = T(k)e ikn for n > 1; 
= A(k) for n = 0; 
= e tkn + R{k)e~ lkn for n < -1. 



(28) 



Here < k < ir. The first line of eq (27) then yields the 
eigenvalue 



i i 



■ cos k. 



2 2 

The next three lines yield the scattering coefficients 

i sin k 



(29) 



A(k) = 7— 
v ' (l-e)e lfc +e-cosfc 

T(k) = (1 - e)A(fc); 
R(k) = (1 - e)A(k) - 1. 



(30) 
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There are also scattering solutions to eq (27) correspond- 
ing to the fictitious particle coming in from the right 



)k _ e -ik 



R(k)e lkn for n > 1; 
= A(k) for n = 0; 

= T(k)e- lkn for n < -1. (31) 

By symmetry the scattering coefficients for this state are 
also given by eq (30). 

There are no bound state solutions to eq (27). The 
scattering solutions we have found all have real eigen- 
values. In principle, since H is non-Hermitian, com- 
plex eigenvalues are also possible. However it turns out 
there are no solutions with complex eigenvalue that are 
biorthonormalisable. It will be seen that the scattering 
solutions we have found constitute a complete set. 

The next step is to find the left eigenvectors that obey 

+ \V r+1 + 1^-1 = *>* \r\ > 2; 



1 + e ,A 



(32) 



Eq (32) is the transpose of eq (27). The left eigenvectors 



are 



^+) k = T(k)e tkn for n > 1; 
= A(k) for n = 0; 

= e ikn + K(k)e- tkn for n < -1. (33) 



and 



^-) fc = e~ ikn + K(ky kn for n > 1; 
= A(k) for n — 0; 



T{k)e 



ikn 



n < -1. 



(34) 



The scattering coefficients are given by 

i(l-e)sinfc 
[> (e-cosfc) + (l-e)e ifc ' 
T(k)=A(k); 
TZ(k) =A(k) - 1. 



(35) 



Having found the left and right eigenvectors, by anal- 
ogy with eq (25), we now posit the completeness relation 

[ g (</4 +)fe *4 +)fc + ^- )fc *^T )fc ) = S mn . (36) 

The proof of this completeness relation, an important 
element of the analysis, is carried out in Appendix A. 

The expansion of the initial correlation vector indi- 
cated schematically in eq (22) may now be written 



Cm(0)= i ^ (+ \k)^ )k +^-\k)4>^ k ]. (37) 
The correlation vector at depth t is now 

c m (t) = £ §A(fc)' [aW(fc)4 +)fc + {k)^ )k 



(38) 



as previously shown schematically in eq (23). 

The expansion coefficients a(k), obtained using the 
completeness relation (eq 36), are 



(+)fc*. 

n ) 



a(+>(fc)= ]T c «(°)^ 

n— — oo 
+oc 



(39) 



as previously indicated schematically in eq (26). To en- 
sure convergence of the sums in eq (39) we set c m (0) — > 
e -ri\m\ and takc ^ o at tne en d. Using cqs (33), (34) 

and (35) we perform the sums exactly to obtain 



oW(t) = oH(k) 



(-)r 



A(k)* + 2A{k)* 



p -ik-T) 



t g—ik — rj I J gi/c — r/ 



ik—rj 



CC 



= 27ryl(fc)*5(A:) + [1 - .4(fc)*]icot 



(40) 



The last line of eq (40) is obtained by taking the limit 

7]^ 0. 

Substituting eq (40) in eq (38) and making use of eqs 
(28), (29), (30), (31) and (35) we finally obtain 



co(t) 



i e r 

1 - e 7T 7 



cos 2 (^>( fc /2) 
e 2 -(26-1) cos 2 ft" 1 J 



Eq (41) is the exact expression for the evolution of Co(t) 
that we sought. 

Finally we would like to re-express eq (41) in terms of 
standard special functions. Some of the manipulations 
will prove useful later in the analysis of injection. Intro- 
duce the z-transform 



11 

c o( z ) = z ' c oW = Y~e~ 
t=o 



.if 

n Jo 



dk 



cos 2 (fc/2) 



1 



e 2 - (2e- l)cos 2 (fc/2) 1 - zcos 2 (fc/2)' 



(42) 
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The fc-integral may now be performed 2 to yield 



CoO) 



1 



(S 2 w(t)) » 4 - - / ds S - 1/2 e-* ln(1+s) (s + S 2 )- 1 
° k Jo 



(1 -£)(!-*) (l-e)(2e-l)(l-a(e)z) 



dss^^e^is + d 2 )- 1 . (47) 



(1 - 2e)vT^z(l - a(e)z) 



(43) 



For brevity a(e) = e 2 /(2e — 1). Upon inversion of the 
z-transform (details relegated to Appendix B) we obtain 



1 1 

5 n Jo 

In the first line of eq (47) we have changed the integra- 
tion variable from x to s — x — 1. Again changing the 
integration variable from s to s = s/6 2 we obtain 

-s 2 tS 2 ' 



(S 2 w(t)) = 



1 



co(t) 



1 1 



1 — e 7T£ 



1 



ds- 



(48) 



(44) 



'Comparing eq (7) we conclude that close to the critical 
point and in the large depth limit (S 2 w(t)} does indeed 
have a scaling form with exponents 



Comparing an integral representation for the hypergeo- 
metric function [E7] 



= 1, y = 2 



F(a, b, c; s) 



r(c) 



T(c-b)T(b) 



dx(x - l) c -"- V" c (x - s) 
(45) 



and scaling function 

T{u) = 1 - 



ds- 



1 



(49) 



(50) 



„2 W = i r(i/2)r(t + 3/2) 

l-e tt£ r(t + 2) 

2e- 1 



valid for Re c > Re b > and |s| < 1, we conclude 

xF\l,t + 3/2,t + 2;=^^\ (46) 

Eq (46) is the final result of this section. It is an exact 
formula for the evolution of fluctuations with depth, in 
terms of known special functions. As a practical matter 
eqs (41) and (44) are equivalent to eq (46) and will prove 
more useful. 

C. Scaling Limit 

Eq (46) gives the exact evolution of load fluctuations 
for the g-model without injection. It is valid for all t and 
all distributions, P(f). From our point of view however 
it is more interesting to examine the scaling limit of large 
depth behaviour near the critical point. 

To derive the scaling limit we start with eq (44) — eq 
(41) would have served just as well — and consider the 
limit t > 1 and 6 = 1 — e — > 0. We do not make any as- 
sumption about the relationship between t and 1/8. We 
obtain 



Fig 4 shows a plot of !F{u). As anticipated the asymp- 
totic behaviour of the scaling function is 



T(u) 



1 



1 



for u — > oo 



>7TU 

for u 



0. 



(51) 



We conclude that the saturation depth scale £ corr ~ 1/S 2 . 
For very great depths t S> £ CO rr> the fluctuations saturate 
to the value 1 /S as found earlier by analysis of the steady 
state (eq 19). For small depths, l<i< £ corr they grow 
as 



(S 2 w(t)) 



2 r 
-r=Vt. 

\/7T 



(52) 



This behaviour must persist at all depths right at the 
critical point as will be explicitly confirmed in section 
III. 

In summary, we have shown that the singular distribu- 
tion is an isolated critical point in the space of g-models. 
There is a (vertical) correlation length that diverges as 
the critical point is approached. We have determined 
the exponents 9 and ip and the scaling function J-(x) 
introduced in eqs (5), (6) and (7). In context of river 
networks we have found that any (/-model with stream 
splitting (hence non-zero 5) has a (possibly very long) 
correlation length in the direction of flow. Such a net- 
work is therefore not scale invariant on sufficiently long 
scales 3 . 



2 Extend the range from to 2ir. Then use the standard trick 
for turning an angular integral into a contour integral around 
the unit circle in the £ plane via the substitution £ — > e lk . 
See for example p 409. 

3 Strictly, to analyse a river network the appropriate initial 
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III. CRITICAL POINT DISTRIBUTION 

Right at the critical point in 1 + 1 dimensions it is 
possible to analyse the dynamics of the entire distribu- 
tion U(w, t). Since there is no vertical length scale at the 
critical point we expect that in the large depth, scaling 
limit 

U(w,t) = ^H(wt T ). (53) 

Eq (53) implies that at the critical point the variance 
should grow as £ _3T_W in the scaling limit t 1. From 
eq (7) we had surmised that the variance would grow as 
f or S = 0. Hence the exponents 0, <p of the preceding 
section and u>, T of this section are not independent; they 
satisfy the relation 3T + oj + 9/ip = 0. Below we calculate 
the exponents u> and T, explicitly verify the exponent re- 
lationship and obtain the scaling function Ti(s). 

Again as an aside to experts we note that the expo- 
nents lo and T might more naturally have been written 
lj — > r/vza, T — > — \jvza. These expressions follow 
from the 5 — > limit of the more general scaling hy- 
pothesis for the entire distribution close to the critical 
point (eq 174). However for this section we have elected 
to make the more restricted hypothesis, eq (53), and to 
give the exponents single character names. 

In this section too we neglect injection. At the critical 
point injection is a relevant perturbation. The form we 
derive is therefore a transient that will break down at 
sufficient depth. Provided the injection is weak however 
that depth could be very great. 

Majumdar and Sire |^8| have analysed the scaling limit 
of Tl(w, t) when injection is present; however it does not 
appear straightforward to take the injection — > limit 
in their expression. It would also be desirable for the 
case of non-zero injection to have a simple explicit for- 
mula for the crossover of n(w, t) from the transient we 
derive (eq 53) to the injection dominated, large depth 
limit. Presumably this can be accomplished by extract- 
ing the suitable limit of the results of ref |2£| , or by direct 
calculation, but we do not attempt it here. 

A. Disorder Average 

As in section II we assume a uniform load is applied to 
the top layer (eq 2). To obtain the distribution U(w,t) 
following ref [|| we consider the quantities 



Z r (p,t) = (ex-pip 5>»(i)> (54) 

ri=l 

where r — 1,2,3,... By translational invariance 
Zx(p,t) = {expipwx(t}) 

OO 

= Y / e vw n(w,t). (55) 

tu=0 

Note that for the critical g-model without injection the 
load on a site is an integer. Thus Zi(p,t) is the discrete 
Fourier or z-transform of the distribution Tl(w, i); p is the 
transform domain variable conjugate to w. Z2(p, t) simi- 
larly encodes the joint probability distribution of load on 
neighbouring sites and so on. 

For the business at hand the imaginary part of Z r (p, t), 

Z r (p,t)=lmZ r (p,t), (56) 

is especially valuable. It is evident from eq (55) that 

OO 

Z 1 (p,t)=^siix(pw)IL(w,t). (57) 

By using Fourier's identity 
2 f n 

— / dk sin kn sin km = 8 mn (58) 
Wo 

and eq (57) we can extract the distribution Tl(w,t) from 
Z\ (p, t) via 

2 r n 

Il(w,t) = - dpsm(pw)Z 1 (p,t) (59) 
k Jo 

for w = 1,2,3,... We cannot obtain H(w — 0, t) in this 
way from Z\ (p, t) , but we can obtain it from the normal- 
isation of H(w, f); 

OO 

U(w^0,t) = l-^U(w,t). (60) 

w — 1 

The benefit of considering the quantities Z r (p, t) is that 
they obey a simple linear evolution equation. Following 
ref B write 

r 

Z r (p, t + 1) = (exp ip w n (t + 1)) 

71=1 

r 

= (exp ip{wi (t)fi (t) +} j w n (t) 
+w r+1 (t)[l-f r+1 (t)}}). (61) 



condition is to load a fraction of randomly chosen sites in 
the t = layer, rather than the uniform load analysed here. 
However we do not expect our conclusion regarding correla- 
tion lengths is sensitive to initial conditions. 
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To obtain the second line we have used the g-model evo- 
lution eq (1). Since the weights in layer t depend only 
on fractions in the preceding layers we can perform the 
average over f\(t) and f r +i(t) separately in eq (61): 

(expiptui (*)/!(*)>/! = i [l + expiptoi(i)] ; 

(expipw r+1 (t)[l - / r+ i])/ r+1 = ^ [1 + cxp ipw r+1 (t)} . (62) 

Substituting eq (62) in eq (61) we obtain the evolution 
equation 

Z r {p,t + l) = ^Z r _ 1 (p,t) + l -z r ( P ,t) + ^Z r+1 {p,t) 

(63) 

where we have again made use of horizontal translational 
invariance. 

Note that eq (63) is linear. Hence it is obeyed sepa- 
rately by the real and imaginary parts of Z. Z therefore 
evolves according to 

Z r (p,t + 1) = iz r _i(p,t) + ^Z r (p,t) + jZ r+1 (p,t) 

(64) 

Eq (64) is reminiscent of a tight-binding lattice 
Schrodinger equation for a free particle on a half-line 
(since the site index r > 1). 

The main results of this subsection are eqs (57) and 
(59) that define the relationship between U(w, t) and 
Z\ (p, t) and eq (64) that controls the evolution of Z r (p, t) 
with depth. 



B. Solution and Scaling Limit 

We wish to solve eq (64) subject to the initial condition 

Z r (p, t — > 0) = sin pr. (65) 

This follows from the assumed uniform load applied to 
the top layer and eqs (54) and (56). Schematically, eq 
(64) has the form 



Z r (p,t + 1) = ^H rs Z s (p,t). 



(66) 



It is easy to verify that our initial condition is an eigen- 
function of H; 



^ H rs sin ps = Q + ^ cos pj sin pr. 



Hence eq (64) has the remarkably simple solution 
Z r (p,t) = Q + ^cospj sinpr. 



(67) 



(68) 



Substituting eq (68) in eq (59) we obtain the desired ex- 
pression for 

2 r /ll \* 

H(w,t) = — / dpsm(pw) - + - cosp I sinp (69) 
n Jo V 2 2 J 

for w = 1, 2, 3, . . . 

The integral over p can be performed exactly by a stan- 
dard contour integration trick (see footnote 2) to yield 



U(w,t) 



1 



(2t) 



4* (t + l-w)\(t-l + w)\ 
1 (2t)\ 



4* (t- 1 -w)\(t + l + w)\ 
for w = 1, 2, . . . , t — 1 

= i m 

4* (t + l-w)\{t-l + w)\ 
= for w > t + 1. 



for w = t, t + 1 
(70) 



We now use eq (60) and (70) to obtain Tl(w -» 0, t). The 
sum proves tractable and yields 



U(w -> 0,t) = 1 



1 (2t+l)l 
4* (t+ 1)!*!' 



(71) 



Eq (70) and (71) are the exact expressions for H(w, t) for 
the critical g-model without injection. 

Much more interesting than the exact formula is the 
scaling limit of large depth. We now assume f > 1 but 
we will make no assumptions about the relative size of w 
and t. To derive this limit we return to eq (69) and write 



1 1 

2 + 2 C ° Sp 



-*p 2 /4 



(72) 



justified (inside the integral) for large t. Hence we obtain 
a Gaussian integral 



U(w,t) = - dpsm(pw)pe- tp2/4 

K J-tt 



4 w 



(73) 



Comparing eq (53) and (73) we see that at large depth 
n has the anticipated scaling form with exponents 



u = 1, T 
and scaling function 



His) = -^se- s . 



(74) 



(75) 



Eq (73) holds for w > 1. In the same large depth limit 

(76) 



II(«;->0,t)«l--?=i=. 

V 71 " \ t 
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The distribution of load thus consists of a spike at zero 
load followed by smooth behaviour for non-zero load 
given by eq (73). At great depths it is extremely prob- 
able that the load on a given bead is zero; most of the 
weight of the distribution is in the spike. 

From the distribution of load, eq (73), it is easy to con- 
firm that its variance (eq 4) grows without bound as the 
square root of depth, as we had earlier inferred from the 
scaling function T (cf eq 52). 

It is instructive that the exact formula, eqs (70) and 
(71), is so cumbersome; the scaling limit, eqs (53), (74) 
and (75), emerges only when we plumb the depths. 

IV. EFFECT OF INJECTION 

In this section we consider the q- model in 1+1 dimen- 
sions taking into account injection. We will assume that 
the weights of the beads are independent and identically 
distributed with mean (7) and variance (SI 2 ). To probe 
the behaviour of the model we will assume that a uniform 
load is applied to the top layer (eq 2). We will study how 
the mean square load (w 2 (t)} evolves with depth since the 
mean load has the trivial variation 

(w(t)) = 1 + (I)t. (77) 

Near the critical point we expect that the mean square 
load should have a scaling form 

(w 2 (t)) = ±C(t5^(SI 2 )5-^(I)5- K )- (78) 

We can guess all the exponents and obtain some fea- 
tures of the scaling function from simple arguments. The 
load on a particular bead at depth t is a random linear 
combination of the weights of the beads in the layer above 
plus a term, due to the applied load, that does not de- 
pend on the weights, /„. Hence the scaling function has 
to be of the form 

{tB » (t)> = i_r m + @£mw) + ^Imn 

(79) 

In the limit of zero injection eq (79) should reduce to 
our result in section II. Thus 

= 1, tp = 2 (80) 

and T has the same form (eq 50) as in section II justifying 
the recycling of these particular symbols. 

By rewriting the average weight at depth t (eq 77) as 
1 + tS 2 (I) /5 2 we conjecture 

k = 2. (81) 



To obtain /i we imagine that the system is very close to 
the critical point. Then for times that are not too long, 
effectively, it will behave as it would right at the critical 
point. At that point the weight of each bead zig-zags 
down lines that merge but do not split. If we add the 
squares of the loads on all the beads on layer t we will 
obtain the sum, over all the beads above layer t, of their 
squared deviation from the average weight (I) plus other 
terms. Hence (%2 n w n (t) 2 ) = (5I 2 )Nt+ other pieces that 
do not depend on (SI 2 ). Here N is the number of beads 
in a layer. By translational invariance we conclude 

(w 2 (t)) w (SI 2 )t+ others. (82) 

In eq (82) "others" represents contributions to (w 2 (t)) 
that do not depend on (SI 2 ). Comparing eq (82) and 
(79) we see that for small values of its argument 

M(u) S3 u (83) 

and the exponent 

M = 1, (84) 

needed to cancel the S dependence at small depths 

With the exponents in hand we can analyse the be- 
haviour of (w 2 (t)) at small depths (compared to 1/S 2 ). 
This behaviour would persist out to all depths right at 
the critical point. For the term independent of injection 
we have already obtained the exact result, eq (52). For 
the term that depends on (SI 2 ) we have just worked out 
the behaviour in this limit, including the precise coeffi- 
cient (eq 82). For the term that is proportional to (I) 
we argue that for small u, K(u) ~ u z / 2 to cancel the 5 
dependence, leading to 

(w 2 (t)) ~ (I)t z ' 2 + others. (85) 

Similarly the contribution of the term that is propor- 
tional to (I) 2 is 

(w 2 {t)) ~ (I) 2 t 5 / 2 + others. (86) 

The last result has a simple interpretation. We have 
seen in section II that without injection at the critical 
point the mean weight at depth t is 1; the mean square 
weight <~ \ft. With injection the average weight at suf- 
ficient depth is ss (I)t. If we assume that uniform in- 
jection does not change the distribution, only its scale, 
then since the mean is inflated by a factor (I)t, the mean 
square should be inflated by a factor (I) 2 t 2 , leading to 
eq (86). The same interpretation can be used to derive 
the behaviour of the last term in eq (79) in the limit 
t ^> 1/S 2 , the opposite of the limit we have so far con- 
sidered. In that limit, in the absence of injection, the 
fluctuations saturate at the value 1/S. Hence we expect 
this term to behave as 

<™ 2 (*)> ~ \(I?t 2 . (87) 
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We can check some of these deductions by making con- 
tact with Majumdar and Sire, who have analysed the 
entire distribution of load at the critical point Q . Fol- 
lowing these authors let us imagine that the injection 
term is very small, with the squared mean (I) 2 , signif- 
icantly smaller than the variance (SI 2 ). According to 
our analysis ultimately the fluctuations at the critical 
point should grow as t 5 / 2 , but the depth at which the 
t 5 / 2 term (eq 86) overtakes the term linear in t (eq 82) 
could be very great; it diverges as l/(/) 4 / 3 . Majumdar 
and Sire arrived at the same value 4/3 for this crossover 
exponent. Moreover, since they argued that right at the 
critical point (the only case they considered) there is only 
one independent exponent, we have made contact with 
their entire analysis as regards exponents. 

In summary we anticipate that near the critical point 
the mean square load will follow the scaling form (eq 
78). Using simple arguments we have conjectured val- 
ues for all the exponents [eqs (80), (81) and (84)] and 
guessed some features of the scaling function. As a check 
we have made contact with the critical point analysis of 
Majumdar and Sire and recovered the known value of the 
crossover exponent, 4/3 p8| . In the remainder of this sec- 
tion we will fully confirm the deductions we have made 
above. We will obtain an exact formula for the evolution 
of the mean square load; the exponents, 9, tp, [i and n; 
and the scaling function, C. 



A. Disorder Average and Exact Solution 

As in section II our strategy is to analyse the evolu- 
tion of the correlation function, c n (t); the mean-squared 
weight (w 2 (t)) — co(t). The analysis is given in outline 
since most of the needed technical elements have already 
been described in section II. Here we shall focus on the 
new complications introduced by consideration of injec- 
tion. 

Following the method of section IIB we first obtain 
the evolution equation for the correlation function, now 
including injection. Schematically this equation has the 
form 



H mn is the same matrix as in eq (15). The effect of injec- 
tion appears in the inhomogeneous term £ m . Explicitly 



U = 2(1) + (2t + l)(I) 2 + (SI 2 )S m=0 . 



(89) 



Our strategy to solve eq (88) is to first expand c(t) and 
£ (t) in terms of the right eigenvectors of H: 



A 



A 



urn" 



(90) 



As discussed before, the expansion amplitudes a\ and £\ 
are calculated by use of the left eigenvectors 



«A 



In section IIB we have calculated the amplitudes for 
c m (t -> 0). We found aW(fc,t-> 0) = d~\k,t -> 0) = 
a(k, t — > 0) with 

a(k, t -> 0) = 27rA(k)*S(k) + [1 - A(k)*]i cot £ (92) 

Here Aik) is given by eq (35). Similarly (k, t — ► 0) = 
£H (k, t 0) = f (Jfc, t -> 0) with 

£(M) = (5I 2 )A(ky + {2(1) + (2t + 1)(I) 2 } 

{2TrA(k)*S(k) + [1 - A(k)*]i cot |}. (93) 

Substituting the expansions eq (91) into the evolution 
eq (88) shows that the dynamics of the amplitudes for 
different right eigenvectors is decoupled and is given by 

a x (t + 1) = Xa x (t) + &(t). (94) 

To solve this dynamics we introduce the z-transforms 



a\{z) = y^QAft); 
t=o 

oo 



t=0 



to obtain 



a\(t -> 0) z£ x {z) 
a\{z) = — ; h 



1 - zA 1 - z\ 
Combining eq (90) and (96) we conclude 

a x (t^0) z£ x {z)~ 



c m (z) = 

A 

More explicitly 

c (z) = 2 



1 - zX 



1 



(95) 



(96) 



(97) 



dk 



10 

+2z 



a(k,t -> Q)A(k) 
1 - zA(fc) 



dk- 

o 1 - zX(k) 



(98) 



Here cq(z) is the z-transform of co(t); A(k) is given by 
eq (30); A(fc), by eq (29); and a(k,t -> 0), by eq (92). 
£,(k,z) is to be obtained by z-transforming eq (93). 

Now all the pieces have been assembled. It remains 
to perform the k integral and invert the z-transform. 
The fc-integrals may be performed by the standard con- 
tour integration method mentioned in footnote 2. The z- 
transforms can all be inverted as illustrated in Appendix 
B. 
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After much calculation we find 
(w 2 (t)) = F(t, e) + (6I 2 )M(t, e) + (I)K(t, e) + (I) 2 L(t, e) 

(99) 

with 

e 1 



F(t,e) 
M(t,e) 

L(t,e) 



1 — e 7re 

-7T^2 + - < ^^nFi + F 2 ]; 
(1 — e) 1 7r 



2e 2 



1-e (1 - e) 3 7re 
e 4 + 2e 3 - e 2 2e 2 



— r[tFi + F 2 ]; 



(1-6)5 (l- e )3 

--^-r[(4t 2 - t)Fi + (8t - 5)F 2 + 8F 3 ]. (100) 

07T£ 



We have put an overline on F(i, e) to avoid confusion with 
a hypergeometric function. For brevity we have written 



r = 



r(i/2)r(t + 3/2) 



r(t + 2) 

3 

2' 



3 2e — 1 

F„ = F ( n,t+-,t + 2; — — 



(101) 



in eq (100). 

Eq (100) is the final result of this subsection. It gives 
the evolution of load fluctuations for the g-model with in- 
jection in 1+1 dimensions. It holds for any distribution 
of fractions P(f) and at any depth. 



B. Scaling Limit 

More interesting than the exact results is the scaling 
behaviour that emerges for t ^> 1 and 6 = 1 — e ^ 0. To 
derive this behaviour it is useful to express the hyperge- 
ometric functions in eq (100) via the integral representa- 
tion, eq (45). The asymptotic behaviour of TFl has been 
analysed in section IIC [cf. eqs (47) and (48)]. The cor- 
responding analysis of rF 2 and rF 3 is very similar and 
finally leads to 

F(t,e)-,i|l-i* 1 («) 

M(t,e) — ^\~l + -u^ 1 (u) + -^ 2 (u)\; 

K{t, e) 1 f 2« + 2 - —$i(u) - -* 3 (u) 

_H [ u 2 $!(u) + 2u$ 2 (u) + 2$ 3 (u)] }. (102) 
Here u = t5 2 . For brevity we have written 



<&„(u) = / cis- — — . 

V ' Jo (l + s 2 ) n 



(103) 



Comparing eq (79) to (102) we conclude that the expo- 
nents are 9 = l,ip = 2, k=1 and fj, = 1 as conjectured. 
It is also straightforward to extract the scaling functions 
!F(u), A4(u),K.(u) and C(u) from eq (102). The scaling 
functions are plotted in figs 4, 5, 6 and 7 respectively. 

The asymptotics of the integrals are analysed 

in Appendix C. Using those results we conclude that for 
small u 



1/2. 



M[u) ps u; 
JC{u) 

C{u) 



8 



30F 
16 
150F 1 



,3/2. 



,5/2 



(104) 



For large w 









- 2 A" 




£(u) s 




r(u) P 


iu 2 . 



(105) 

Substituting the small u asymptotics in eq (79) we ob- 
tain the behaviour for depths small compared to 1/5 2 . 

(W 2 (t)) « A t V2 + {Jja>t + {J) 8 t 3/2 + (/)^ t5 /2. 
V 71 " 3vtt 15V7T 

(106) 

This behaviour would persist for all depths right at the 
critical point. Note that eq (106) agrees with the forms 
conjectured in eq (82), (85) and (86) (including the nu- 
merical coefficient in the first case) . It is hardly necessary 
to add that eq (106) is consistent with the critical point 
analysis of Majumdar and Sire since it leads, by the ar- 
guments given earlier in this section, to their crossover 
exponent 4/3 |§]. 

The large u asymptotics give the behaviour at depths 
large compared to 1/8 2 . We find 

(w 2 (t)) = 1 + (SI 2 )~ s ^=t^ + (J)±2i + (I) 2 -/. (107) 

The term proportional to (I) 2 has the form anticipated 
in eq (87); at the greatest depths this term is dominant. 

In summary we have shown that the singular distribu- 
tion is an isolated critical point. Near the critical point 
the fluctuations in load have the scaling form eq (79). 
We have derived this scaling form and all the exponents. 
The results are in agreement with expectations based on 
simpler (non-rigorous) arguments. 
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V. HIGHER DIMENSIONS 



w(n, t + 1) = ^ fft(n — u, t)w(n — u,t). 



(109) 



We now turn to the g-model in D + 1 dimensions. The 
quantum Hall multilayer and river networks are both 1+1 
dimensional systems; bead-packs however are described 
by the 2+1 dimensional g-model. The behaviour of the 
model as a function of D is of intrinsic interest moreover. 
We will find that right at the critical point the growth ex- 
ponents vary smoothly with dimension for D < 2. Above 
D = 2 they become fixed, revealing D = 2 as the upper 
critical dimension for the critical case. Off the critical 
point we expect the fluctuations to grow according to 
a scaling function T(x) (eq 7). We will study how the 
function and exponents vary with dimensionality below 
D = 2. For simplicity in this section we neglect injection. 



Eq (109) is the D + 1 dimensional generalisation of cq 
(!)• 

The fractions for a particular bead are assumed to be 
drawn from a distribution that is symmetric with respect 
to direction and respects the constraint eq (108). It fol- 
lows 



We write 



(fu) ~ 2D' 



</!> = 



i 



+ 



22D 2 2D 



(110) 



(111) 



A. Model and Disorder Average 

First we must generalise the description of the g-model, 
so far confined to 1+1 dimensions. The case of 2+1 di- 
mensions is easy to visualise. Fig 8 illustrates a square 
lattice composed of two interpenetrating square sublat- 
tices. The co-ordinates of sites n — (m, 712) are both even 
for the black sites; both odd for the grey. The displace- 
ments from a site on either sublattice to its four nearest 
neighbours on the other sublattice are (±1, ±1). We will 
denote these displacements u. In the g-model planes of 
such square lattices are stacked vertically. The beads al- 
ternately occupy only even or odd sublattices. Denoting 
the depth of a layer t, for t even only the even sublattice 
is occupied; for t odd, only the odd sublattice. Viewed 
in three dimensions the beads occupy a body-centered 
cubic stucture. In the same sense, Fig 1 can be viewed 
as a body-centered square structure. 

Now consider a D dimensional simple cubic lat- 
tice. The co-ordinates of a site are specified by n = 
(m, n 2 , ...,«£)) where rii are integers. For the even sub- 
lattice the 7ij are even; for the odd sublattice, they are 
odd. Each site has 2 D nearest neighbours on the other 
sublattice. We denote the displacements (±1, ±1, . . . , ±) 
to these neighbours u. The D + 1 dimensional g-modcl 
consists of D dimensional cubic lattices stacked in the 
"vertical" t direction. In alternate t slices only the even 
or odd sublattices are occupied by beads. 

It is assumed that a random fraction of the load on 
each bead is supported by its 2 D neighbours in the layer 
below. The fractions must sum to one; 



5> 



i. 



(108) 



where e is a parameter that characterises the distribution 
of fractions. From the sum constraint eq (108) it follows 



(fu 1 fu 2 ) — 



1 



1 



2 2D (2 D -1)2 2D ' 



(112) 



for Mi + 1 u 2 . The fractions are assumed to be indepen- 
dently and identically distributed for different beads. 

For the singular distribution all the fractions are zero 
except one. The probability for each fraction to be one 
is 1/2 D . It is easy to calculate e = 2 D — 1 for the sin- 
gular distribution using eq (111) and to verify eq (110) 
and (112) are satisfied. Since e = 2 D — 1 for the singular 
distribution we shall use 5, defined by 

(113) 



(2 D - 1)' 

as our measure of the distance of a distribution from the 
critical point. 

As before it is useful to consider the correlation func- 
tion 



c(m, t) = V^(u)(n, t)w(n + m, t)). 



(114) 



Note that m is a D dimensional vector with even inte- 
ger entries for both t even and t odd. The correlation 
function therefore lives on a simple cubic lattice in D di- 
mensions. By rescaling, as in section IIA, we reduce the 
lattice constant of this lattice to one so that the compo- 
nents of the vector m are now integers. The variance in 
load is related to the on-site correlation by 



(5w 2 (t)) =c(m-»0,t)-l. 



(115) 



Here /g is the fraction of load transmitted by the bead 
to the neighbour separated by a horizontal displacement 
of u. Hence the dynamics of the g-model is governed by 



Following the discussion of section IIA and using eqs 
(109), (110), (111) and (112) it is easy to show that the 
correlation function evolves with depth according to 
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c(m, t+l) = ^cim, t) + c (™ + *>> f ) 



2 D + 1 



7+2 E c ( 



For example, c(p, t — > 0) denotes the Fourier transform 
of c(rh,t) at t = 0; no z-transform is implied. 

Performing both transforms on eq (116) we obtain 



2 D+2 



z{m + b,t) + 



c(p, z) = c(p, t->0) + zc(p, z)S{p) 



-^d E C (™ + M) + ■^ 5 c{m,t)8rh=o 



+ 



(2 D - 1) 



(2^-1) 



2 D 



2 D - 1 



zc(to — > 0, z). 



(120) 



^ = i{i+~E e ^+i E e ^ 



6-nnn 



6— nn 

6-nnn 1 * ' ' 1 / j ^ J 

J ^ b— n...n 

i)^D E C (^ + M)<W6=0 ( U6 ) _ (1 + COSpi) (l + COS P2 ) (l+COS PD ) 



(2^ - 1) 



6— n. . .n 



While reading eq (116) it is useful to recall that the cor- 
relation function lives on a Z?-dimensional cubic lattice. 
For D = 2 each site has four nearest neighbours and four 
next-nearest neighbours. For general D, each site has 2D 
nearest neighbours; 2 2 C(D,2) next nearest neighbours; 
2 3 C(L>, 3) third nearest neighbours; and 2 D C(D 1 D) D th 
nearest neighbours. In eq (116) b denotes the displace- 
ment from a site to any of these neighbours; nn denotes 
nearest neighbour; nnn, next nearest; and so forth. 

In the next subsection we will solve eq (116) for c(m — > 
0, t) subject to the initial condition that a uniform load 
has been applied to the top layer. Thus c(m, t — ► 0) = 1 
for all to. 



B. Solution 



It is easy to verify that 

c(m, t — > oo) = — for to = 0, 
o 

= 1 otherwise, 



(117) 



is a steady state solution to eq (116). Eq (117) shows 
that the variance (5w 2 ) saturates at sufficient depth in 
all dimensions for all distributions except the singular. 

We now calculate the evolution of the variance with 
depth using a method different from that of section II 
p9[ . First we z-transform the (discrete) t dependence of 
the correlation function, 



c(m, z) 



t=o 



c(m, t)z , 



(118) 



and Fourier transform the space dependence, 

c(p, z)=J2 e~^" l c(m, z). (119) 

m 

The use of the same symbol for the correlation and its 
transforms, although customary, is potentially confusing. 



(121) 



is a "structure factor" for the cubic lattice. It will also 
prove convenient to define 



G(p, z) = 



1 



1 - zS{p) ' 



(122) 



Both S(p) and G(p, z) have helpful physical interpreta- 
tions that we shall make use of below. For the moment 
we rearrange eq (120) to obtain 

c(p, z) = c(p, t -> 0)G(p, z) 

+ (1 - 5)zc{m -> 0,z)[l - S(p)]G{p,z). (123) 

By inverting the Fourier transform we can turn eq (123) 
into an expression for c(m — ► 0,z). After further re- 
arrangement 



c(to — > 0, z) 



l-(l-S)zJ J ^[l-S(p)}G(p, z y 



(124) 



Eq (124) is a general expression for c(to — * 0, z) for an 
arbitrary initial condition. For uniform loading of the 
top layer 



c(p,i-> 0) = (2ir) D 8(p). 



(125) 



It follows from eq (121) and (122) that G(p -»■ 0, z) 
1/(1 — z); hence eq (124) simplifies to 



c(m — ► 0, z) = (1 — z)~ 



x<l-(l-8)z 



dp 



(27T) 



D 



\l-S(p)]G(p,z) 



(126) 



Eq (126), together with the definitions of the structure 
factor (eq 121) and G(p, z) (eq 122), constitutes an exact 
formal evaluation of the variance with depth. To obtain 
(8w 2 (t)} explicitly it only remains to peform the integral 
over p and to invert the z-transform. We return to this 
task in the next subsection. We conclude this subsection 
with a useful interpretation of S(p) and G(p, z). 



14 



Eq (116) with e — > resembles the Schrodinger equa- 
tion for a particle on a .D-dimcnsional cubic lattice with 
hopping to the nearest neighbours, the next nearest 
neighbours, and so on to the D th nearest neighbours. It is 
not difficult to sec that the cigcnstates of this Schrodinger 
equation are plane waves. S(p) is the dispersion relation, 
the eigenvalue at wave vector p. From eq (121) we see 
that the energy level spectrum is a continuous band be- 
tween zero and one. 

The momentum space Green's function for this tight- 
binding lattice would normally be written 



G(?,E) = 



1 



E-S(p)- 



(127) 



Comparing eq (127) to eq (122) we see that G(p, z) is es- 
sentially the Green's function with E — > 1/ z. It is famil- 
iar from quantum mechanics that the real space Green's 
function at the origin, 



g{m^ 



■ E) =J 



dp 



1 



{2ir) D E-S{p)' 



(128) 



regarded as a function of (complex) E, has a branch cut 
running from E — to E — 1 , the interval that supports 
the eigenvalue band. It is not difficult to use the famil- 
iar arguments to conclude that, regarded as a function of 
complex z, c(m — > 0, z) has a branch cut along the line 
z = 1 to oo (onto which the segment [0,1] maps under 
the transformation E — > 1/z). The analytic properties of 
c(m — > 0, z) will prove useful in the next sub-section. 



C. Scaling Limit 

In this subsection we study the evolution of the vari- 
ance in the large depth scaling limit. Thus t ^> 1 and S 
is zero or very close to it throughout. 

An advantage of studying the large depth limit is that 
we do not have to calculate c(m — > 0, z) exactly; it is only 
necessary to calculate the leading behaviour as z — * 1. 
One way to understand this is to consider the critical 
case <5 = 0. In this case we expect that at great depth 



c(m -> 0, t) ~ t x 



(129) 



It is easy to show that for f(t) = t x , the z-transform 
is T(x + 1)/(1 — z) x+1 plus less singular terms. Thus 
for a function that behaves as t x for large t also the z- 
transform is 



T(x + l) 

(1 - z y+ l 



less singular. 



(130) 



If we know the leading singularity of c(m — > 0, z) as z — > 1 
we can use eq (130) to read off the large depth behaviour. 

Another way to see that we only need the behaviour 
of c(m — > 0, z) as z — > 1 is to consider inverting the z- 
transform by the contour integral method of Appendix 



B. This is accomplished by folding the contour over the 
branch point of c(m — > 0, z) at z = 1 and integrating 
along the cut. In that integral c(m — > 0, z) is weighted by 
a factor that decays extremely rapidly away from z = 1 
at large depths. 

Our goal therefore is to analyse the z — > 1 behaviour 

of 



G(*) 



dp 



1 



(2tt)^ 1 - zS(p) 



(131) 



since by a straightforward re-arrangement the integral in 
eq (126) simplifies to 

/ jSr [l SmG{P ' Z) = { l -~z) G{Z) + \ (132) 

Insight into the behaviour of G(z) can be gained by ex- 
panding S(p) around p = to obtain 



G(z) 



dp 



1 



(2tt) d (l-z)+p 2 /4' 



(133) 



If we set z = 1 in cq (133) the integrand diverges as 
p — > for D < 2; it is regular in more than two dimen- 
sions. Thus in more than two dimensions G(z) has a 
branch point at z = 1 but there is no actual divergence. 
In two dimensions or less there is an actual divergence. 

The leading behaviour of G(z) above two dimensions 
is thus simply obtained by setting z = 1 in cq (131): 



G(z) w G(l) for D > 2. 



(134) 



In two dimensions we can obtain the singularity by recog- 
nising G(z) to be a Jacobi elliptic integral. Square lattice 
Green's functions are known to be related to Jacobi's 
elliptic functions; but since our lattice features next- 
nearest neighbour hopping, in addition to the custom- 
ary nearest neighbour hopping, we outline the analysis 
in Appendix D. The result is that for z — > 1 



G(z) ss --ln(l - z) for D 

IT 



2. 



(135) 



For D < 2 we obtain the singular behaviour of G(z) in 
Appendix D. The result is 



G{z) 



r(i-D/2) 



(l-z) 



D/2-1 



for D < 2. 



(136) 



An important feature revealed by this calculation is that 
the singular behaviour of G(z) is controlled by the long 
wavelength behaviour of G(p, z) for all D < 2; it breaks 
down as D — > 2. Although it is instructive to do the cal- 
culation for continuous D to examine the D — > 2 limit, 
the only case that is physically relevant is of course the 
integer dimension D = 1. 
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Equipped with the leading behaviour of G(z) in all 
dimensions we now obtain the long time behaviour of 
(Sw 2 (t)). At the critical point we set S = and substi- 
tute eqs (132), (134), (135) and (136) in eq (126). Except 
in two dimensions the z-transforms may be inverted by 
inspection of eq (130). For two dimensions we must re- 
sort to the method of Appendix B and finally obtain 

(Sw 2 (t)) = n / 2 - 1 sin (jpj t D ' 2 for D < 2 

= — for D = 2 
In* 

= ^ for D > 2. (137) 

As indicated by the simple steady state solution, at the 
critical point the fluctuations grow without bound as a 
power of t for all dimensions. The exponent becomes in- 
dependent of D for D > 2 revealing D = 2 as the upper 
critical dimension. 

By substituting eq (132) and (136) in eq (126) we can 
also obtain the behaviour of (5w 2 (t)) away from the crit- 
ical point for less than two dimensions. Inverting the 
z-transform by the method of Appendix B we find 

(Sw 2 (t)) = ^HtSn- (138) 
Here the exponents 

9 = 1, p = -| (139) 
and the scaling function 

1 2 f°° 1 us 2/D 

^ = d-^dL " s TT^ cxp -^ (140) 

with qn — T(l — D/2) sin(7rD/2)/Y / 7r D a dimension de- 
pendent constant. Again, only the result for D = 1 is 
physically meaningful; in this case eq (140) coincides with 
the result of section II. 

In summary the main results of this section are that for 
all distributions, except the singular, at sufficient depth 
the load fluctuations saturate and (in agreement with ex- 
periment) there are no horizontal correlations in load (eq 
117). The saturation value of the load variance diverges 
as the critical point is approached. At the critical point 
the load fluctuations grow without bound as a power of 
depth (eq 137). Below two dimensions this exponent de- 
pends on dimensionality; above two dimensions it is con- 
stant, revealing D = 2 as the critical dimension. At the 
critical dimension the growth of fluctuations is tempered 
by a logarithmic factor as might be expected at a critical 
dimension. We have also evaluated the scaling function 
that describes the growth and saturation of load fluctu- 
ations near the critical point for D < 2. 



VI. QUANTUM HALL MULTILAYER 
A. Models 

In this section we turn to the chiral wave models that 
are believed to adequately describe the surface electronic 
states of a quantum Hall multilayer. We begin by exam- 
ining the circumstances under which the quantum net- 
work model of Saul, Kardar and Read Q discussed in 
section I is equivalent to a g-model. 

Following Saul, Kardar and Read, the first step is to 
identify pairs of links (joined by vertical grey bars in Fig 
9) as "beads" . The "load" on a bead is the total proba- 
bility that the electron is on either of its two constituent 
links. Load propagates from left to right now rather than 
top to bottom as it did in our earlier depictions of the 
g-model. For this reason we will label the vertical co- 
ordinate n and the horizontal co-ordinate t here (see fig 
3). 

To analyse how load propagates consider an elemen- 
tary vertex of the Saul, Kardar and Read model shown 
in Fig 9. The wave function amplitudes are related via 

(*)-»(*)■ 

here S is a random 2x2 su(2) rotation matrix. Saul, Kar- 
dar and Read assumed the S'-matrices were drawn from 
the invariant distribution for the su(2) group f3(i|| . The 
loads on beads A, B and C are respectively l^il 2 + 1^2 1 2 , 
|(/.i| 2 + |0 2 | 2 and|03| 2 + ^4| 2 • Byunitarity \^i\ 2 + \^ 2 \ 2 = 
1 4>3 1 2 + 104 1 2 - Thus bead A sends a fraction / of its load 
to neighbour B and the remainder 1 — / to neighbour C. 

A key feature of the Saul, Kardar and Read model is 
that the distribution of the fractions, P(f) is indepen- 
dent of the input amplitudes tpi and ip2- This follows 
from the assumed group invariant distribution for the S'- 
matrices. It is this feature that allows the Saul, Kardar 
and Read model to be mapped onto a g-model. 

To derive the distribution of the fractions recall that 
an su(2) matrix may be parametrized S — xq + ix.a 
with (xq,x) real and subject to x^ + x 2 — 1. If we take 
ipi = 1, ip2 = then / = Xq + x\. From the invariant 
distribution for su(2) matrices, 

P(x ,x) = -S{x 2 +x 2 -l), (142) 

it is not difficult to show that the fraction / follows the 
uniform distribution, P(f) = 1 for < / < 1. 

Now suppose the wave function is known through the 
vertical slice t = 0. We could propagate the wave func- 
tion t slices to the right using the quantum Saul, Kardar 
and Read model. Alternatively we could calculate the 
load in the initial layer and propagate it to the right us- 
ing the q-model with uniform distribution. Either way 
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the load we obtain in layer t would be the same statis- 
tically. This is the sense in which the Saul, Kardar and 
Read model is equivalent to the g-model. 

Note that the g-model does not keep track of phase 
information. The mapping is useful only under circum- 
stances that the phase information is unimportant. Be- 
low wc will discuss some problems of wave packet dy- 
namics for which the mapping is useful. The mapping 
can also be used to study vertical transport in the quan- 
tum Hall multilayer in the limit of large circumference 
but we do not discuss that application here. 

An obvious circumstance when the phase information 
is important and the mapping cannot be used is if peri- 
odic boundary conditions are imposed in the horizontal 
i-direction, as would be appropriate for a multilayer in 
the fully phase-coherent, mesoscopic regime. Phase in- 
formation is needed to match the wavefunction after it is 
propagated around the circumference. We will develop 
this point in a more technical way in subsection VI C. 

Another case in which a quantum network model will 
map onto a g-model is if the wavefunctions and S- 
matrices are chosen to be real and the S'-matrices are fur- 
ther assumed to be distributed over the subgroup of ro- 
tations about the y-axis with appropriate invariant mea- 
sure. The fraction distribution P(f) = (l/7r)/ _1 / 2 (l — 
J)- 1 / 2 f or the q— model that results. For most distribu- 
tions of the S'-matrix however it is not possible to obtain 
even the limited mapping between the quantum network 
model and the classical g-model obtainable in this and in 
the Saul, Kardar and Read case. 

Finally we present a convenient continuum model of 
the multilayer surface governed by the Schrodinger equa- 
tion 

- 4iW = +m;_ l (t)i_i(i). (143) 

at 

Since the equation is first order in t, given the wave- 
function at a fixed t slice we can use it to propagate the 
wavefunction to the right, just as in the discrete network 
model. In the transverse direction the model is discrete 
and second -order. Disorder is incorporated by taking the 
hopping elements m n (t) to be random. For a discussion 
of the relationship between onsite and hopping disorder 
see ref j24y26[] - Evidently this model cannot be reduced 
to a classical g-model. 



B. Wave- packet dynamics 

In this section wc briefly discuss wave packet dynamics 
for the models of the previous section. Mathematically 
this problem is identical to the motion of a wave-packet 
in a crystal with noise (temporal randomness). It also 
bears formal resemblance to the directed polymer model, 
an important problem in statistical mechanics. Hence 
it is a problem of general interest and has been studied 



since at least the 1970s from various points of view (see 
ref [|| and refs therein). A considerable amount is now 
known. 

For the Saul, Kardar and Read model wave-packet dy- 
namics can be studied using the mapping to the q-model; 
indeed the mapping was introduced for this purpose. In 
this section we will formulate the problem and summarize 
known results. These results reveal that the g-model and 
the continuum wave model introduced in the last section 
behave in qualitatively similar ways. 

Consider an electron localized at n — at t = 0. This 
wavepacket can be propagated to the right using eq (143). 
As it propagates it will broaden and its mean position 
will deflect. It is interesting to know how the breadth 
and deflection grow with displacement and to analyse the 
distribution of "load" at sufficiently great displacement 
that a steady state is reached. 

The root mean square width of the wave-packet grows 
as the square root of the displacement. This was derived 



for the continuum model in the 1970s |31 1 and it is easy to 
show that the same form is obtained in the Saul, Kardar 
and Read model. The root mean square deflection grows 
as the fourth root of the displacement. This result has 
been obtained numerically and analytically for both the 
Saul, Kardar and Read [6 29 3^] and continuum models 

IIS 

To compare the distribution of load, for the continuum 
model we define the load on an edge as w n (t) = \ip n (t)\ 2 . 
The asymptotic distribution of load, Tl(w, t oo) was 
obtained by Coppersmith et al. for the q- model 0. For 
various distributions of the fractions, P(f), they found 
that n(u>, t — > oo) decayed exponentially with w with 
a power law prefactor that depended on the distribu- 
tion P(f). For the uniform distribution the prefactor 
was a constant. The corresponding result for the contin- 
uum wave model was obtained by ref 25 by mapping 
the problem onto an su(l,l) quantum ferromagnet. Here 
too the result for the load distribution is an exponential 
with a prefactor linear in w. 



C. Field Theory Formulation 

We have emphasized above that the equivalence be- 
tween the Saul, Kardar and Read model and the q-model 
is useful only when open boundary conditions are im- 
posed in the horizontal i-direction; it breaks down for pe- 
riodic boundary conditions needed to describe transport 
in phase-coherent multilayers. The importance of bound- 
ary conditions is also reflected in field theory formulations 
of these models. In ref p^j the continuum model with 
open boundary conditions was mapped onto a Heisen- 
berg ferromagnet. In contrast, with periodic boundary 
conditions a mapping to a supersymmetric analogue of 
the Heisenberg ferromagnet was obtained in refs [E3 26 
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In this section we derive the supersymmetric spin rep- 
resentation following the operator methods of ref 
This derivation highlights the role of boundary condi- 
tions, the feature we wish to emphasize here. It only 
makes use of operator methods and is in this sense more 
elementary than the functional methods of ref [^6| . More- 
over mappings to supersymmetric spin models have re- 
cently been used fruitfully not only to study the multi- 
layer but also to provide non-perturbative insights into 
various other problems of electron localization p4|-|36|| . It 
is hoped that the present derivation, with its emphasis 
on boundary conditions 4 and its use of operator methods 
will prove of interest in this broader context also. 

1. Fermion Representation 

We wish to evaluate G(n, t; n' , t'), the Green's function 
for the continuum model governed by the Schrodinger 
equation, 

d 

- i—G(n, t; ri, t') = m„ (t)G(n + l,t; ri, t') 
ot 

+m* n _ 1 (t)G{n- n',t') 
-i5(t-t')S nn <, (144) 

and subject to the periodic boundary condition 

G(n, t + T; ri, t') = G{n, t; ri, t'). (145) 

Here T is the period in the ^-direction. In ref pi] , 
the Green's function was calculated subject to the chiral 
boundary condition, G(n,t;n' ,t') = for t < t' , leading 
to a simpler field theory formulation. 

The key idea is to reinterpret the co-ordinate t as time. 
Eq (144) then describes a particle on a one-dimensional 
lattice with noise. In second quantised notation the 
(time-dependent) Hamiltonian that governs the motion 
of this fictitious particle is 

ff«(t) = £ [m„(i)cf C * +1 + m^ity^c^} . (146) 

n 

Here c^' creates a Fermion at site n; annhilates it. 
The reasons for the superscript on the Fermion Hamilto- 
nian and on the creation and annhilation operators will 
become apparent shortly. 

The S'-matrix for this model obeys 

-i±S$(t) = H*(t)S§(t) (147) 



subject to S$(t -» 0) = 1. From S%~ l S§ = 1 it is easy 
to verify the useful result 

i^S«- 1 =S«- 1 (t)H F <(t). (148) 

We define 

c*(t) = S*- 1 (t)c n S*(t) (149) 

and similarly for c^(t). 

Now by analogy with finite temperature field theory 
we write the Green's function 

G(n,t;n',t') = Tr [S$(T)<X{t)c$tf)]/Z*(T) 
for t > t' 
= -Tr [S*(T)c$(ir)c*(t)]/Z*(T) 
for t < t'; 

Zp(T) = Tr [Sp(T)]. (150) 

Zp(T) is analogous to the partition function in finite 
temperature field theory. It is easy to verify that G obeys 
the differential eq (144) by making use of eqs (147) and 
(148) and the commutation relation 

= -m„(t)c« +1 - m* (151) 

However 

G(n,T;n',t') = Tr [Sp(T)Sp^ 1 (T)c^Sp(T)c^} (t')]/Zp(T) 
= Tr[S*{T)c B J{t') C K]/zK{T) 
= -G(n,0;ri,t'). (152) 

Thus G obeys antiperiodic rather than periodic bound- 
ary conditions. This problem is fixed by adding a term 
to the Hamiltonian 

ff£(t)^2#(t) + (153) 

n 

Alternatively we may replace Tr — ► STr in eq (150). By 
STr we mean the trace of an operator over all states with 
an even number of fermions minus the trace over states 
with an odd number of fermions. 

We also need an expression for the complex conjugate 
of the Green's function since our ultimate purpose is to 
calculate the disorder average of \G(n, t; n' , t')\ 2 , the dif- 
fuson propagator. To this end we complex conjugate eq 



4 For the effect of boundary conditions on the supersymmetry 
mapping for models such as the Chalker model of the quan- 
tum Hall transition see refs and J3fj where the random 
hopping model in one-dimension is analysed with periodic and 
open (scattering) boundary conditions respectively. 
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(144) to obtain the differential equation obeyed by G* . 
Comparison to eq (144) reveals that we should consider 
A fermions governed by the Hamiltonian 



£ [ m ™(*)< 



At r A 



+ m*(t)tf<£ . (154) 



G*(n, t; n', t') is then given by the right hand side of eq 
(150) if we replace R -> A and Tr -> STr. 

As might be expected the Hamiltonian for the A 
fermions is related to that for the R fermions via a par- 
ticle hole transformation. This symmetry between the R 
fermions and the A holes leads to an su(2) symmetry in 
the fermion sector of the complete field theory formula- 
tion that we obtain below (eq 171). It is also at the root 
of the supersymmetry of the field theory formulation. 

In summary the Green's function with periodic bound- 
ary conditions may be generated from the second quan- 
tized Hamiltonian, H${t) [eq (153) and (146)] using the 
definition eq (150). The complex conjugate of the Green's 
function may be obtained similarly using the Hamilto- 
nian H$(t) (eq 154). Eq (150) and its A fermion ana- 
logue provide exact formal expressions for the Green's 
function for a particular realization of the random tun- 
neling m n (t). These expressions are not particularly con- 
venient to average since m n (t) appears in both numerator 
and denominator. 



n 

The main result of this subsubsection is eq (156). 
It provides a formal bosonic expression for the exact 
Green's function for a particular realization of random 
tunneling, m n (t). A similar expression for G* may be 
obtained by working with the Hamiltonian eq (157). Like 
their fcrmionic counterparts these bosonic expressions 
are not particularly well suited for averaging over dis- 
order. 



3. Supersymmetry 

We now develop an expression for the diffuson suitable 
for averaging over disorder. In Appendix E it is shown 
that 



Z*(T)Z§(T) = 1. 



(158) 



Thus we consider a model that includes both A and R 
fermions and bosons governed by the Hamiltonian 

#susyW - H${t) + H${t) + H§(t) + Hi(t) 

= J2[m n (t)A n +m* n (t)Al]. (159) 



2. Boson Representation 

Alternatively we could interpret eq (144) as a time de- 
pendent Schrodingcr equation for bosonic particles on 
a one dimensional lattice. The corresponding "time"- 
dependent bosonic Hamiltonian in second quantized no- 
tation is 

H%(t) = ]T KWW+i + <-i(*) W-i] • (155) 

n 

Here creates an R boson at site n; b„ annhilates it. 
The Green's function is now defined as 

G(n,t;n',t') = Tr [S*{T)b R n {t)b*} {t')]/ Z*{T) 
for t > t' 
= Tr [Si(T)b^(t')b^(t)]/Z§(T) 
for t < t'; 

Z§(T) = Tr S§(T). (156) 

Here S§ is the bosonic S-matrix and Zg(T) is the bosonic 
analogue of the partition function. 

For greater rigour we must regulate the traces to en- 
sure convergence but for brevity we do not discuss this 
explicitly here. 

The complex conjugate of the Green's function is gen- 
erated similarly if instead of the R bosons we consider A 
bosons governed by 



Here 

Ji n — C n C n+1 C n+1 C n -f D„ 0„ + i O n+1 O n . 

The corresponding S-matrix obeys 

- *^Ssusy(*) = H SVSY (t)S S u SY (t) (161) 

subject to <SsuSY(i — > 0) = 1. A formal solution to eq 
(161) is given by 

SsusyW = Pcxp (i f dhHsvsYih)) . (162) 



Here P is the chronological ordering operator. 
Hence the diffuson is given by 

|G(n,t;n',0| 2 = STr [S av MT)c*(t)c*(t)c$ 
for t > t' 

= STr [S SVSY {T)4ht')c*ht')${t)c£(t)] 
for t < t' (163) 

The content of eq (163) is that to calculate the diffuson 
we must create or annhilate a pair of R and A fermions 
(depending on the time order). Then we must propa- 
gate this state in accordance with -ffsuSY and perform 
an S-matrix weighted trace. The Hamiltonian -ffsusY is 
non-interacting but it is random and time dependent. 
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Eq (163) is an exact formal expression for the diffuson. 
Note the lack of a denominator, eliminated by virtue of 
eq (158). This feature allows us to perform the average 
over disorder easily. For example, 



(SsusyW) = exp[-H S usYt] (164) 



with 



n 

Here we have assumed that the tunneling m n (t) is a 
Gaussian white noise process with zero mean and vari- 
ance 



{m^(t)m^>(t')) = DS(t - t')5 nn ,5 ap . 



,08) 



(166) 



Here rrin\t) = real part of m n (t): rrin\t) — imaginary 
part of m n (t). 

Recall that for a single Gaussian random variable y, the 
phase average (e w ) = e~^ y Eq (164) is analogous to 
this result but with the added complications that Ssusy 
is an ordered exponential, not a simple exponential, and 
the average is over a random process rather than a single 
random variable. To derive eq (165) it is simplest to ex- 
pand the time ordered exponential (eq 162) and average 
term by term. 

Proceeding in this manner we obtain an expression for 
the average diffuson 

(\G(n, t; n', t')\ 2 ) = STr {exp[-H slJS Y (T - t + t')]c R c A 



x exp[-7Y SUSY (i 
for t > t'. 



(167) 



A similar expression may be written for the case t < t'. 
The content of eq (167) is that to calculate the aver- 
age diffuson we must create (or for the other time order, 
annhilate) a pair of R and A fermions and propagate 
the resulting state according to the effective Hamiltonian 
7isuSY- In contrast to Hsusy the effective Hamiltonian 
is not time dependent or random but it is interacting. 

This completes our formulation of the continuum di- 
rected wave model of section VI A as a superspin field 
theory. The main results are the superspin Hamiltonian 
(eq 165) and eq (167) which shows how interesting corre- 
lation functions are calculated in the superspin formula- 
tion. The usefulness of this formulation depends on the 
extent to which the superspin model can be analysed. 

In the remainder of this section we discuss the form and 
symmetry of the superspin Hamiltonian (eq 165). To this 
end it is helpful to introduce special notation for the bo- 
son and fermion bilinears of which Hsusy is composed. 
We denote the fermion bilinears 



J_ 



c A c R 



J x 

J X 



iJy 



J z = ^c R + c A h A -l) 

J = L c m c R_ c At c A + 1) 



(168) 



the boson bilinears, 



K + = b R ^b A ^ = K x + iK y , 
K_ = b A b R =K X - iK y , 

K z = ^b R + b A ^b A + l), 
K=-{b R ^b R -b A ^b A -l); 



(169) 



and the mixed bilinears, 



Mi = ft fl t c R ) m 2 = b A h A , 



(170) 



In eqs (168), (169) and (170) the site indices have been 
suppressed for brevity. In terms of these bilinears we may 
write 

^susy = — 2£>^ ( J n+ ±.J n + J n+ \J n — J„J 

n 

(Kn+i-Kn + K n+ iK n + K n } 

n 

+^EKM 1) + M >i 2) + hc ) 

n 

+^E( £ «ll4 1) +4T 1 4 2) +hc). an) 

n 

Here he denotes Hermitian conjugate and K n+ i.K n — 



It is instructive to study the commutation relations for 
bilinears at the same site n (bilinears at different sites 
simply commute or anticommute). It is easy to verify 
that J+, J_ and J z satisfy angular momentum or su(2) 
commutation relations and J commutes with the other 
three. Similarly K +1 K~ and K z satisfy the su(l,l) or hy- 
perbolic angular momentum algebra — essentially the an- 
gular momentum algebra but with a sign change for the 

K commutes with the other 



K + ,K- commutator 38 



three. The anticommutators of Li, Lj, Mi and Mj are lin- 
ear combinations of the K's and J's. The commutators 
of the J's or K's with the L's or M's are linear combina- 
tions of the L's and M's. Hence these bilinears constitute 
a superalgebra. The J's and K's are commuting elements 
of the superalgebra; the L's and M's, anticommuting el- 
ements. The superalgebra is called u( 1,1 — 2). It includes 
the Lie algebras su(2) and su(l,l) as subalgebras. 

Further insight into the superalgebra is obtained by 
considering the Hilbert space at a single site. This is a 
direct product of the four dimensional fermion space and 
the infinite dimensional two-boson space. The fermion 
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space may be decomposed into irreducible representa- 
tions of the su(2) algebra. The fermion vacuum and the 
state with both R and A fermions present constitute a 
doublet or spin 1/2 representation; the two states with 
one fermion present are singlets. The boson space simi- 
larly decomposes into an infinity of infinite dimensional 
irreducible representations of the su(l,l) algebra 5 . The 
single site Hilbert space thus decomposes rather simply 
into irreducible representations of the direct sum of the 
su(2) and su(l,l) algebra. These subspaces do not con- 
stitute a representation of the whole superalgebra. The 
anticommuting elements mix different irreducible repre- 
sentations of su(2) and su(l,l). In particular they mix 
representations with different spins — a celebrated feature 
of supersymmetry. It is not difficult to decompose the 
single site Hilbert space into blocks irreducible under the 
superalgebra; however this would carry us too far afield. 
More details on the superalgebra are given in ref (26) and 
refs therein. 

Finally we define 

>7tot=X^«- (172) 

n 

Here J denotes any element of the superalgebra such as 
J +1 K Zl Li etc. After some algebra we find 

[WsusY,Jtot] = (173) 

revealing the supersymmetry of the field theory formula- 
tion. 



VII. SUMMARY AND CONCLUSION 

Much of this paper is concerned with the behaviour 
of the g-model close to the critical point. To probe this 
behaviour we imagine that a uniform load is applied to 
the top layer. As the load propagates downward fluctu- 
ations develop in the distribution of load. Coppersmith 
et al. [Q studied the entire distribution of load at very 
great depth where it was presumed that a steady state 
had been reached. In contrast we study only the variance 
of the distribution of load but we analyse its evolution 
with depth. Our purpose is to study this evolution for 
all distributions of the fractions, P(f), particularly those 
close to the singular distribution (the critical point). 



In section II we consider the (/-model in 1+1 dimen- 
sions without injection (the weight of the beads is ne- 
glected). In this case the average load does not vary with 
depth since the total load is the same in every layer; it is 
merely redistributed by the q- model dynamics. For the 
growth of the variance, by analogy to critical phenomena, 
we make the following hypotheses: For all distributions 
P(f) except the singular distribution we posit that the 
variance will saturate at sufficient depth. Both the satu- 
ration depth and the saturated variance are expected to 
diverge as the distribution approaches the singular dis- 
tribution. We introduce S, a measure of the distance of 
a distribution P(f) from the singular distribution, and 
conjecture that the saturation depth £ corr will diverge as 
1/(5 V ; the saturated variance, as l/5 e . More specifically, 
we expect that close to the critical point the variance 
will have a scaling form, eq (7). For the singular dis- 
tribution we expect that the variance will grow indefi- 
nitely as a power of the depth. Close to the critical point 
and at depths shallow compared to the saturation depth 
the variance should grow as it would right at the critical 
point. From this and from eq (7) we deduce a relationship 
between the critical exponents 8 and tp and the exponent 
that describes the growth of the variance right at the 
critical point; namely we expect that at the critical point 
the variance will grow as t 8 /^ . In sections IIB and IIC 
we derive an exact formula for the variance as a function 
of depth (eq 46) and study its scaling limit (t ^> 1, S — ► 
but with t5 v arbitrary). These calculations bear out all 
the expectations enumerated above, provide the precise 
form of the scaling function [eq (50) and Fig 4] and yield 
the exact exponents (eq 49). 

In section III we characterise the critical point more 
fully by analysing the evolution with depth of the entire 
distribution of load right at the critical point in 1+1 di- 
mensions. In the absence of injection the critical point 
is a simple model of random walkers that coalesce upon 
contact; hence it is quite straightforward to derive these 
results. We present them because they illuminate the re- 
sults of the previous section. At large depth it is found 
that the distribution of load consists of a large spike at 
zero load together with a smooth part [eq (73) and (76)]. 
It is overwhelmingly probable that the load on a bead 
is zero; most of the weight of the distribution is in the 
spike. The smooth part follows the anticipated scaling 
form (eq 53). Its width grows as the square root of the 



5 Let \n + m, n) denote a state with (n + m) R-bosons and 
n A-bosons on the site. The infinite dimensional subspace 
with m a fixed integer and n = 0,1,2,..., for m > 0, or 
n — —m, —m + 1, = m + 2, . . ., for m < 0, is invariant under 
the four K operators. These subspaces corresponding to dif- 
ferent values of m constitute the irreducible representations 
of the su(l,l) algebra. 
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depth, consistent with the exponent found in section II to 
describe the growth of the variance of load at the critical 
point. 

In section IV the effect of injection is included. For 
simplicity we consider only 1+1 dimensions. We assume 
that the weights of the beads are independent and iden- 
tically distributed random variables. The behaviour of 
the mean load is still not very interesting. It grows lin- 
early with depth (eq 77). Close to the critical point we 
conjecture that the variance will have the form eq (78). 
We are able to deduce all the exponents in eq (78) and 
to obtain some limiting behaviours of the scaling func- 
tion through simple (non-rigourous) arguments. These 
conjectures are all verified by the exact calculation of 
sections IV A and IV B which provides the precise form 
of the scaling function [eqs (99), (102) and (103)] and 
yields all the exponents [eqs (80), (81) and (84)]. We find 
that beyond a crossover depth the variance (normalized 
by the squared mean) saturates. The saturation value 
and the crossover depth both diverge as the critical point 
is approached. At depths less than the crossover depth 
the variance grows as it would right at the critical point 
(eq 106). The behaviour at the critical point has many 
crossovers if the weight of the beads is small compared to 
the applied load. In this case at first the variance grows 
as the square root of the depth as it was found to do in 
section II in the absence of injection. At greater depths 
there are crossovers to growth as t and t 5 / 2 , as first the 
effects of large rare fluctuations in the weight of a bead 
and then mean injection assert themselves. Ultimately at 
the critical point the variance grows with the 5/2 expo- 
nent but the depth at which this behaviour sets in can be 
very great if the mean injection is small. This depth di- 
verges as (J) -4 / 3 . The crossover exponent 4/3, deduced 
by simple arguments and then via exact calculation in 
section IV, agrees with the value previously obtained by 
a different method by Majumdar and Sire ]^j|. In their 
work Majumdar and Sire only study the behaviour right 
at the critical point. However at this point they calculate 
the dynamics of the entire distribution of load whereas 
we study only the variance. 

In section V we turn to the g-model in D+l dimen- 
sions. For simplicity we neglect injection in this section. 
We find that right at the critical point the variance grows 
as a power of depth in all dimensions except two (eq 137). 
The power is given by D/2 for D < 2. For all dimensions 
above two the growth is linear. This shows that D = 2 is 
the upper critical dimension for this problem. For D = 2 
we find a linear growth of the variance tempered by a log 
factor as might be expected at the critical dimension. 

An intriguing feature of the critical behaviour we ob- 
tain is that it is exhibited at all. For ordinary continuous 
phase transitions the renormalization group provides a 
framework to understand the critical behaviour. We are 
not aware of any such framework for the q-model. 

Random critical points are notoriously difficult to anal- 



yse in general. The feature that allows us to analyse 
the g-model is that the two point load correlation func- 
tion [defined by eqs (8) and (114)] evolves with depth 
according to a simple linear equation. In section II we 
analyse the evolution by expanding in the eigenvectors 
of an appropriate linear operator. There are some sub- 
tleties posed by the non-Hermiticity of the linear oper- 
ator, making it necessary to prove that its eigenvectors 
are complete (further complicated by the infinite dimen- 
sionality of the vector space). Nonetheless we like this 
approach because it parallels transfer matrix methods 
used for equilibrium critical phenomena. We find that 
the large depth scaling behaviour is controlled by the 
low energy long wavelength eigenfunctions of the non- 
Hermitian "Hamiltonian" . Another virtue of this ap- 
proach is that with about the same effort it yields both 
the variance and the correlation functions. However we 
have left analysis of the correlation functions open for 
later work. Here we focus entirely on the variance of the 
load. In section V we analyse the variance using another 
technique based on transform methods. 

Our analysis, neglecting injection, confirms that the q- 
model has essentially no horizontal correlations in the 
steady state for any distribution except the singular. 
This agrees with experiments on bead packs. The bulk 
of our results however are concerned with the q-model 
close to the critical point. Bead pack experiments such 



as those of ref |12| appear to be far from the critical 
point. We estimate 5 0.5 for this experiment. It is 
not obvious how to tune the parameter for bead packs 
to access the critical behaviour we analyse here. Claudin 
et al. have also studied the horizontal steady state cor- 
relations of the q-model without injection away from the 
critical point H . They employ a continuum limit and ar- 
rive at conclusions similar to eq (19) in section IIB. The 
main focus of their work however is to explore a tensor 
model of stress propagation in granular matter, intended 
to supplant the q-model. 

Interpreted in terms of river networks our results show 
that allowing a small amount of river splitting in a Schei- 
degger network introduces a length scale in the vertical 
direction. On sufficiently long length scales such a net- 
work is not scale invariant. This resembles the finding of 
Narayan and Fisher |2fl] . In their model too there was a 
parameter that controlled river splitting. Their networks 
were not scale invariant unless river splitting was tuned 
to zero. However their model appears to be in a differ- 
ent universality class as its vertical correlation exponent 
is different from the value ip — 2 we obtain here in sec- 
tion II. Presumably the difference is because their rule 
for stream splitting was non-local and depended on the 
entire history of the network upstream from the split. 

Taken together with the model of Narayan and Fisher 
it appears that river splitting is a perturbation that spoils 
the scale invariance of river networks. It is therefore in- 
teresting to ask whether such networks exist in Nature. 
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River deltas are one possibility. Traced backwards they 
may constitute networks of merging streams that occa- 
sionally split. Even for river basins it might be interest- 
ing to examine the extent to which streams split. In this 
context it is worth noting that some of the data against 
which river scaling laws are tested are based not on ac- 
tual maps of the river network but on networks that are 
indirectly inferred according to certain rules from digital 
elevation data obtained from satellite images. The rules 
by which the network is inferred from the elevation maps 
exclude the possibility of splitting [Q . 

In summary the g-model is rich in applications and 
behaviour and yet analytically tractable by elementary 
means a combination of circumstances that invites fur- 
ther exploration. Among the many problems that remain 
open we conclude by mentioning two: For the q-modcl be- 
yond the saturation depth there is no correlation in the 
horizontal direction but in the vertical direction there 
are very strong and long ranged correlations |59|. We 
have not obtained the precise form of these vertical cor- 
relations for the g-model either in steady state or at the 
critical point. It would be very interesting (and straight- 
forward) to obtain these forms and the crossover between 
them. Second it would be interesting to obtain the dy- 
namics of the entire distribution of load near the critical 
point. We have not attempted to do this except right at 
the critical point. 

A natural scaling hypothesis is that the full distribu- 
tion of load, neglecting injection, will be of the form 

n(w,t,S) = W - T Q(wS 1/a ,tS l/z ). (174) 

The exponents in eq (174) are r = 2, vz = 2 and a = 1. 
Their values are fixed by our result for the variance away 
from the critical point derived in section II and the re- 
sult for the entire distribution at the critical point 5 = 
derived in section III. We also know that for x — ► and 
y —>■ 0, the presently unknown function Q has the asymp- 
totic behaviour 

eo*,tf)« -^4m e ~ x2/y ( 175 ) 

to be consistent with the critical point distribution (eq 
73) derived in section III. 

In the second part of this paper we turn to chiral wave 
models that are believed to describe the surface electronic 
states of a quantum Hall multilayer. In section VI A we 
discuss circumstances under which the quantum network 
model of Saul, Kardar and Read (described in the in- 
troduction) is equivalent to the g-model. In section VI 
B we compare known results about the behaviour of the 
g-model to a continuum chiral wave model that cannot 
be mapped onto a g-model under any circumstance. The 
two are found to behave in qualitatively similar ways. 

A circumstance under which the mapping to the q- 
model is not useful is when periodic boundary conditions 



must be imposed in the chiral direction. Physically this 
is because of the interference of electron paths that wind 
around the quantum Hall multilayer. Such long range 
interference cannot be captured by the classical g-model. 
In this phase-coherent or mesoscopic regime, the chiral 
wave model has been studied via a mapping to a super- 
symmetric spin model [^3j,^6|. In section VI C we derive 
this mapping in a way that emphasizes boundary con- 
ditions. Our derivation makes use of operator methods 
and is hence more elementary than the derivation of rcf 
p6| that makes use of mixed functional integrals over 
Grassman and bosonic variables. We do not attempt 
further analysis of the superspin model here; the inter- 
ested reader should consult papers on multilayer trans- 
port, particularly refs |^(| and |Q that provide a nice 
overview of the early work on this problem. 

Mappings to superspin models have been useful not 
only in the study of the quantum Hall multilayer but 
have also recently lead to new non-perturbative results 
and insights into other important problems of electron 
localization |34| ^] . Hence it is hoped that our deriva- 
tion, with its emphasis on boundary conditions and use 
of elementary operator methods, will be of interest in this 
general context. 

Note added: While writing this paper we learnt of an 
e-print by Rajesh and Majumdar on spatio-temporal cor- 
relations in the Takayasu model and the g-model [ ElJ . 
These authors derive many interesting results comple- 
mentary to ours. In this paper we concentrate on the 
behaviour close to the critical point. For the q-modcl 
Rajesh and Majumdar concentrate on length scales long 
compared to our vertical correlation length, £ C orr; the 
crossovers and scaling functions that we study are tran- 
sients that are invisible in their asymptotic formulae. On 
the other hand they have derived both vertical and hor- 
izontal load correlation functions; this paper is limited 
(in practice but not in principle) to the study of the vari- 
ance of load. Among their interesting findings are (i) 
They find power law correlations in the vertical direction 
both at the critical point and away from it addressing in 
part a question raised above, (ii) They emphasize the in- 
teresting structure of the horizontal correlation function, 
including injection, at great depth. 

Although the goals are a bit different, there are points 
of intersection between the two papers with regard to 
technique. Rajesh and Majumdar too exploit the lin- 
earity of the relation that describes the evolution of the 
correlations with depth and solve it using the method 
of section V. An overlapping result is a formula for the 
variance at the critical point in 1+1 dimensions includ- 
ing injection. At the large depths studied by Rajesh and 
Majumdar the last term in our eq (106) should dominate. 
Rajesh and Majumdar obtain the same exponent 5/2 and 
the same numerical pref actor 16/15-^/7? providing a nice 
check on both calculations. 
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APPENDIX A: PROOF OF COMPLETENESS 

First let us recall the principles of biorthogonal ex- 
pansion (see, for example, ref ^7j, p884). We discuss the 
simplest case of a finite NxN dimensional non-Hermitian 
matrix H mn . Consider its eigenvectors 



E 



(Al) 



In context of biorthogonal expansion these eigenvectors 
are called the right eigenvectors. For simplicity we will 
assume that the right eigenvalues are non-degenerate in 
this case. The bad news regarding the right eigenvec- 
tors is: (i) A may be complex, (ii) There is no guarantee 
that there are N eigenvectors (needed to span the vec- 
tor space), (iii) Eigenvectors corresponding to different 
eigenvalues are not necessarily orthogonal. 

Now consider the left eigenvectors, defined as the 
eigenvectors of W . (i) If A is a right eigenvalue then A* is 
a left eigenvalue (Proof: The coefficients for the charac- 
teristic polynomials of H and H> are complex conjugates 
of one another), (ii) There are as many left eigenvectors 
as right, (iii) Left eigenvectors are orthogonal to right 
eigenvectors. 

The last point merits elaboration. Let ipn denote the 
left eigenvector with left eigenvalue A*. Thus 



n 

According to (iii) above 



(A2) 



(A3) 



Eq (A3) is the biorthogonality relation. It may be proved 
by noting 

mn n 

= A'£>»)V£ (A4) 

n 

whence £„(V£)*^' = for A ^ A'. 

In general there is no guarantee of completeness, but 
in this case assume that N eigenvectors have been found. 
Then we can prove the completeness relation 



A 



(A5) 



The proof follows from the observation that if there are 
N eigenvectors, any vector a n may be expanded as 



a ^ 



(A6) 



Completeness then follows from biorthogonality, eq (A3). 

The problem in section IIB presents some complica- 
tions not present in the pedagogical discussion above. 
Among them are degeneracy, an infinite dimensional vec- 
tor space and a continuous spectrum. Nonetheless the 
broad strategy is the same. In section IIB we found left 
and right eigenvectors and we conjectured biorthogonal- 
ity and completeness relations. To justify the analysis 
of section IIB we must prove the completeness relation. 
That is the purpose of this appendix. Note that we can- 
not simply assume completeness is true — because the ma- 
trix H is non-Hermitian there are no theorems to guar- 
antee it. Nor can we prove completeness by counting 
eigenvectors as in the finite dimensional discussion above. 

The proof of completeness is remarkably simple and 
direct. We substitute the exact expressions for ipm^ k 
and 4>n that we have derived, eqs (28), (31), (33) and 
(34), on the right hand side of eq (36) and verify the 
completeness relation by explicit evaluation of the inte- 
gral. There are nine cases to consider corresponding to 
n = 0, n > 0, n < and m = 0, m < 0, m > 0. 

For illustration we analyse the case of n = 0,m = 0. 
We must evaluate 



- / dkA*(k)A{k) 



(A7) 



where A(k) and A(k) are as given in eqs (30) and (35) 
Since the integrand is symmetric in k we extend the range 
of integration from — n to ir and substitute z 
obtain a contour integral about the unit circle 



e lk to 



dz 1 



(l-e)(z + lf 



2ni z e 2 (z - l) 2 - (1 - e) 2 (z + l) 2 



(A8) 



Evaluation via Cauchy's theorem reveals that the integral 
equals one as required for completeness. 

The remaining eight cases also succumb to this method 
of analysis. 



APPENDIX B: INVERSE Z-TRANSFORM 

Consider the series /(i), t = 0, 1, 2, . . . Its z-transform 
is defined as 



f(z) = J2f(t> 



(Bl) 



t=o 
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Some z-transforms can be inverted by inspection. For 
example the inverse transform of (1 — az)^ 1 is evidently 

(1 - az)- 1 - f(t) = a*. (B2) 

In other cases the inverse transform can be found by per- 
forming the complex integral 

The contour C must enclose the origin but no singulari- 
ties of f(z). 

For illustration let us analyse 

f(z) = (l-z)- 1 ' 2 {l-az)- 1 (B4) 

needed to go from eq (43) to (44) in section IIB. Here 
a > 1. f(z) has a pole at 1/a and a branch cut at z = 1 
(see fig 10). We deform the contour C that encloses the 
origin to contours C\ and Ci that encircle the pole and 
pass above and below the branch cut. Hence obtain 

(B5) 

The first term is the contribution of the pole; the second, 
of the branch cut. 

APPENDIX C: ASYMPTOTICS OF $ N (U) 

The asymptotics of the functions defined by eq 

(103) are needed to obtain the asymptotic behaviour of 
the scaling functions in sections IIC and IVB. 

The large u behaviour poses no difficulty. Evidently 

pK 1 

$ n( M )~"V^ as li^OO (CI) 

for all n. The small u behaviour is a bit more subtle. 
Moreover, it turns out that due to cancellations we will 
need as many as five or six terms in the small u series for 
$„ to obtain the leading behaviour of the scaling func- 
tions. 

For definiteness consider the small u behaviour of 

* l( „)=y o ds^-^. (C2) 
The leading term is obtained by setting u = 0, 

*i(0) = |. (C3) 

To obtain the next term it is tempting to expand the 
integrand in powers of u but this leads to divergent inte- 
grals. The divergence signals that the asymptotic series 
is not a simple power series in u. 



It turns out the next term goes as yfu. To show this, 
and to efficiently obtain many more terms in the series, 
consider 

9( X ) = l dsj^. (C4) 

We will show that g{x) is regular about x = and that 
its asymptotic behaviour is a simple power series. To this 
end we observe that g(x) obeys the first order differential 
equation 

^g - 2xg(x) + ^ = 0. (C5) 

x = is a regular point for this equation; hence we at- 
tempt a series solution 

g(x) = b + hx + b 2 x 2 + ... (C6) 

We find b\ = —\Pr and the simple recurrence relation 

b n = -b n -2- (C7) 

n 

Evidently b = g(0) = n/2. Hence we obtain the asymp- 
totic series 

9{x) = g - V™ + 2 X ~ 3^ x + 4 ~ 15^ X + ' ' ' 

(C8) 

Substituting x — > \fu we conclude 

= | - V^u 1/2 + ^u- ^V^u 3/2 + ^u 2 
~^V^u 5 ' 2 + ... (C9) 
for small u. Similarly 

$ 2 («) = | - | u + ^V^u 3/2 - ^nu 2 + ^Vtvu 5/2 + ... 



APPENDIX D: LATTICE GREEN'S FUNCTION 

1. Two Dimensions 

Consider the Green's function in two dimensions for 
the lattice Schrbdinger equation discussed in section VA. 
The real space Green's function at the origin is given by 

(Dl) 
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[cf. eq (121) and (128)]. We consider real E > 1. In this 
appendix we show that 



(D2) 



Here K is a complete elliptic integral of the first kind. 
From the well-documented properties of these integrals 
or by direct analysis of eq (D8) below it follows that as 
E^ 1+ 



(D3) 



In section V C we are interested in the behaviour of G(z), 
eq (131), as the real variable z — > l - . Comparing eq 
(131) to (Dl) we see that 



G(z) 



-G E 



Hence the singularity of G(z) as z — > 1 is 



G(z) 



— ln(l - z) 

7T 



(D4) 



(D5) 



Eq (D5) is the main result of this section of the Appendix. 

To demonstrate eq (D2) we regard p as a complex vari- 
able p — > x + iy. The integral over p in eq (Dl) may 
be regarded as an integral around the contour sketched 
in Fig 11 since the two vertical segments cancel by the 
periodicity of the integrand and the horizontal segment 
at infinity makes no contribution because the integrand 
vanishes along it. The integrand in eq (Dl) has a simple 
pole at p = iy, where y satisfies 



(D6) 



cosh (T X E 



cos(fc/2)' 



with residue 



iE\ 1 - 



cos z 



E 



Hence by Cauchy's theorem 



(D7) 



yy ' Ej_„2n\ E 



-E C ° S 2 



-1/2 



(D8) 



Comparing to the definition of the elliptic integral of the 
first kind 



K{k) 



we obtain eq (D2). 



tt/2 



d6(l - k 2 sin 2 (9)~ 1/2 



(D9) 



2. Below Two Dimensions 

In this section we analyse the singular behaviour as 
z — ► 1 of G(z) in less than two dimensions. The ap- 
proximate long wavelength expression for G, eq (133), 
provides a useful starting point. 

To analyse the divergence in D = 1 we would note that 
the integrand in (133) is a sharply peaked Lorentzian. 
This justifies working to quadratic order in S(p) and ex- 
tending the range of integration (strictly confined to the 
Brillouin zone, — tt < p < tt in one dimension) to ±oo. 
Result: G{z) = (1 - z)- 1 ' 2 . 

To continue this result to non-integral D we use 't 
Hooft and Veltman's dimensional regularization trick 
El. We write 



G(z) 



ds 



dp 



exp-.s[(l-z)+p 2 ]; (D10) 



extend the range of integration, outside the Brillouin 
zone and over all p-space; and replace 



dp 



n 



D 



since the integrand in eq (D10) is isotropic in p. Here 
£Ie> = 2 v / 7r' D /r(D/2) is the total solid angle in D dimen- 
sions (some familiar special cases: fl\ — 2, = 2n, Q3 = 
4?r, fi 4 = 2?r 2 .) The result is 

T(l - D/2) 



dpp 



D-l 



(Dll) 



G{z) 



-(1 



\D/2-l 



(D12) 



for D < 2. 

This analysis breaks down in two dimensions and 
higher because the integrand diverges as p — » 00. The 
divergence is an artifact of the quadratic approximation 
in eq (133) and of extending the integral outside the Bril- 
louin zone. The spurious divergence is revealed in eq 
(D12) as a pole in the Gamma function factor as D — > 2. 



APPENDIX E: ANALYSIS OF PARTITION 
FUNCTION 

The purpose of this appendix is to show that the par- 
tition functions for bosons and fermions cancel. Thus 



Tr [Sf{T)]Tt [S§(T)] = 1. 



(El) 



A similar relation holds for the advanced bosons and 
fermions. We discuss the retarded case explicitly. For 
brevity the superscript R will be omitted. 
We write the fermion S-matrix as 



S F (t) = exp J2 c n c ^ S fW>- 



(E2) 



iSi?(t) is then governed by the Hamiltonian eq (146) with- 
out the extra term included in eq (153). 
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To make further progress we introduce e™(£), the solu- 
tion to the Schrodinger eq (144) 

- iJUf(t) = mi(t)ef +1 (t) + mU(t)eU(t) (E3) 



subject to ef(t — > 0) = <5„;. 
The scattering formula 



(-1 



(E4) 



will prove very useful. To derive it, rewrite eq (E4) as 

Sp^c^t) = e l n {t) C[ (E5) 
i 

and regard it as an ansatz with the functions e n (t) un- 
specified. Making use of eqs (147), (148) and (151), the 
t derivative of the left hand side is 

S F (ty 1 {m n (t)c n+1 + m*_ 1 (i)c r[ _i}<Sir(i) 

- ci{mi(t)e l n+1 (t) + mU(ty n -i(t)}- (E6) 
i 

To obtain the second line we have made use of the ansatz 
(E5). Comparing eq (E6) to the t derivative of the right 
hand side of eq (E5) we conclude that e l n (t) does obey 
the Schrodinger eq (E3). This completes the proof of the 
scattering formula (E4). 

Another relation that will prove useful is 



S F (t)\0) = |0>. 



(E7) 



This follows because the Hamiltonian (eq 146) annhilates 
the vacuum; S F is the (chronologically ordered) exponen- 
tial of the Hamiltonian. 

Equipped with these results we write the fermion par- 
tition function as 

Z F (T) = Tr { exp [i ir ^ cjcn] S F (T) } 

n 

= (0\S F (T)\0) 



1 



% S F (T)ci\0) 



"1 ,"2 



3j E (VK^niCnsSpiTyi^Cn^O) 
n 1 ,n 2 ,n3 



+ ■ 



(E8) 



The trace is taken over the entire Fock space including 
states with different total numbers of fermions. The al- 
ternating signs are due to the factor exp[i7r^ n c}jC n ] in 
the trace. The factorials are because the sums over the 
site indices rii are unrestricted; hence each state gets 
counted a multiple number of times. 



We now shift the S-matrix to the left using the scat- 
tering formula (E4), make use of the adjoint of (E7) and 
calculate the vacuum expectations of the fermion oper- 
ators (Wick's theorem). The result for the second-order 
term is 



Z F (T) = \ 



-i 2 



-others. 



(E9) 



Fig 12(b) shows a diagrammatic representation of this 
term. Note that the diagram series for the partition func- 
tion Z F (T) contains both connected and unconnected 
graphs. By familiar arguments p7| we can write 



Z F (T) = exp[-fi(T)] 



(E10) 



where the "free energy" f2(T) has the linked diagram ex- 
pansion shown in fig 12(c). 

We turn now to the boson partition function. 

The boson scattering formula 



b n S B (t) =Y J Z l n{t)S B {t)b l 



(Ell) 



can be proved in the same way as eq (E4). Eq (E7) re- 
mains true when we replace Sp — > Sb 

The boson partition function is therefore given by 

Z B {T) = Tr {S B (T)} 
= (0\Sb(T)\0) 

+ Y^Ks B (T)bi\o) 



+ ... (E12) 

This equation resembles eq (E8) but there is an extra 
subtlety in the combinatoric factors. In the two-boson 
case for example, for the offdiagonal terms (ni ^ n-i) the 
factor (1/2!) is to offset double counting as in the fermion 
case. For the diagonal terms, that vanish in the fermion 
case, there is no double counting but the factor (1/2!) is 
needed for normalization. 

By shifting Sb(T) to the left by use of the scattering 
formula we see that the series for Zb(T) is the same as 
for Z F (T) except for the minus signs. Hence 



Z B (T) =exp(+fi(T)) 



(E13) 



where f2(T) is defined by the diagram series in fig 12(c). 

Eq (E10) and E(13) together lead to (El), the result 
we sought to prove here. 
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Fig 1. The q-model of stress propagation through a 
bead pack in 1+1 dimensions. The beads are assumed 
to sit on a regular lattice. Each bead is supported by its 
two nearest neighbours in the layer below. 
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Fig 2. Scheidegger's model: For the singular q- model 
the load zig-zags down lines that merge but do not split. 
In Scheidegger's model of river basins these lines are in- 
terpreted as tributaries merging to form a river. 
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(a) (b) t ► 

Fig 3. A quantum Hall multilayer: Layers of two 
dimensional electron gases are stacked vertically and a 
strong perpendicular magnetic field is applied. The im- 
portant electronic states are at the edge of each layer. 
These chiral edge states propagate in the direction shown 
in (a). A quantum network model for the surface of the 
multilayer is shown in (b). 




Fig 4. The scaling function T{u) describes the growth 
of load fluctuations with depth for the 1+1 dimensional 
(j- model close to the critical point [see eqs (7) and (50)]. 
Injection is neglected. 




Fig 5. Growth of load fluctuations with depth for the 
1+1 dimensional q- model close to the critical point. The 
scaling function M. (u) gives the contribution due to fluc- 
tuations in the weight of beads [eqs (79) and (102)]. Here 
M.(u)/y/u is plotted as a function of u. 



K (u) /u 




Fig 6. Growth of load fluctuations with depth for the 
1+1 dimensional q- model close to the critical point. The 
scaling function JC(u) gives a contribution proportional 
to the average weight of beads [eqs (79) and (102)]. Here 
K(u)/u is plotted as a function of u. 
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Fig 7. Growth of load fluctuations with depth for the 
1+1 dimensional g-model close to the critical point. The 
scaling function C{u) gives a contribution proportional to 
the square of the average weight of beads [eqs (79) and 
(102)]. Here C(u)/u 2 is plotted as a function of u. 




Fig 8. Horizontal slice through the 2+1 dimensional 
g-model. Beads occupy the even (black) or odd (grey) 
sublattice in alternate layers. Each bead is supported by 
its four nearest neighbours in the layer below. 
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Fig 9. Elementary vertex of the Saul, Kardar and Read 
model. 
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Fig 10. Contours for inverting the z transform. 
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Fig 11. Contour for evaluation of two dimensional 
Green's function. 



Fig 12. Feynman diagrams for the partition function, 
(a) Diagram representation of the propagator ef(T). (b) 
Second-order diagrams for the partition function, (c) 
Lowest order diagrams in the infinite series for the free 
energy. 
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